Overview

Setup

(hide)

Click on tabs to display additional information.

Libraries

# Load required libraries for statistical analysis and table creation
library(dplyr)
library(ggplot2)
library(gt)
library(broom)
library(car)
library(emmeans)
library(multcomp)
library(MASS)  # For negative binomial regression
library(rstatix)
library(coin)
library(phyloseq)
library(microViz)
library(tidyverse)

Plotting

# Define treatment order and color palette
treatment_order <- c(
  "A- T- P-",  # Control
  "A- T- P+",  # Parasite
  "A+ T- P-",  # Antibiotics
  "A+ T- P+",  # Antibiotics_Parasite
  "A- T+ P-",  # Temperature
  "A- T+ P+",  # Temperature_Parasite
  "A+ T+ P-",  # Antibiotics_Temperature
  "A+ T+ P+"   # Antibiotics_Temperature_Parasite
)

# Custom color palette matching treatment order
treatment_colors <- c(
  "#1B9E77",  # A- T- P- (Control)
  "#D95F02",  # A- T- P+ (Parasite)
  "#7570B3",  # A+ T- P- (Antibiotics)
  "#E7298A",  # A+ T- P+ (Antibiotics_Parasite)
  "#66A61E",  # A- T+ P- (Temperature)
  "#E6AB02",  # A- T+ P+ (Temperature_Parasite)
  "#A6761D",  # A+ T+ P- (Antibiotics_Temperature)
  "#666666"   # A+ T+ P+ (Antibiotics_Temperature_Parasite)
)

# Create named vector for color scale
treatment_color_scale <- setNames(treatment_colors, treatment_order)

Functions

# Function to extract sample data as dataframe from phyloseq object
samdatAsDataframe <- function(ps) {
  samdat <- phyloseq::sample_data(ps)
  df <- data.frame(samdat, check.names = FALSE, stringsAsFactors = FALSE)
  return(df)
}

# Function to rename variables in phyloseq object
ps_rename <- function(ps, ...) {
  ps <- microViz::ps_get(ps)
  df <- samdatAsDataframe(ps)
  df <- dplyr::rename(.data = df, ...)
  phyloseq::sample_data(ps) <- df
  return(ps)
}

# SourceFolder function
source(here::here("Code", "R", "Functions", "StartFunctions", "sourceFolder.R"))

# Import all helper functions found in `/Functions`
sourceFolder(here::here("Code", "R", "Functions", "StartFunctions"), T)
## 
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
## 
##     discard
## The following object is masked from 'package:readr':
## 
##     col_factor
## The following object is masked from 'package:coin':
## 
##     pvalue
## Loading required package: future
## 
## Attaching package: 'future'
## The following object is masked from 'package:survival':
## 
##     cluster
## 
## Attaching package: 'ggpubr'
## The following object is masked from 'package:microViz':
## 
##     stat_chull
## Loading required package: grid
## ========================================
## ComplexHeatmap version 2.18.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
## 
## If you use it in published research, please cite either one:
## - Gu, Z. Complex Heatmap Visualization. iMeta 2022.
## - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
##     genomic data. Bioinformatics 2016.
## 
## 
## The new InteractiveComplexHeatmap package can directly export static 
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
## 
## This message can be suppressed by:
##   suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
## Loading required package: data.table
## 
## Attaching package: 'data.table'
## The following objects are masked from 'package:lubridate':
## 
##     hour, isoweek, mday, minute, month, quarter, second, wday, week,
##     yday, year
## The following object is masked from 'package:purrr':
## 
##     transpose
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
## Loading required package: doParallel
## Loading required package: foreach
## 
## Attaching package: 'foreach'
## The following objects are masked from 'package:purrr':
## 
##     accumulate, when
## Loading required package: iterators
## Loading required package: GUniFrac
## Registered S3 method overwritten by 'rmutil':
##   method         from
##   print.response httr
## Loading required package: magrittr
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:purrr':
## 
##     set_names
## The following object is masked from 'package:tidyr':
## 
##     extract
## Loading required package: rlang
## 
## Attaching package: 'rlang'
## The following object is masked from 'package:magrittr':
## 
##     set_names
## The following object is masked from 'package:data.table':
## 
##     :=
## The following objects are masked from 'package:purrr':
## 
##     %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
##     flatten_raw, invoke, splice
## Loading required package: vegan
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.6-4
## 
## microbiome R package (microbiome.github.com)
##     
## 
## 
##  Copyright (C) 2011-2022 Leo Lahti, 
##     Sudarshan Shetty et al. <microbiome.github.io>
## 
## Attaching package: 'microbiome'
## The following object is masked from 'package:vegan':
## 
##     diversity
## The following object is masked from 'package:scales':
## 
##     alpha
## The following object is masked from 'package:ggplot2':
## 
##     alpha
## The following object is masked from 'package:base':
## 
##     transform
## Loading required package: ape
## 
## Attaching package: 'ape'
## The following object is masked from 'package:ggpubr':
## 
##     rotate
## The following object is masked from 'package:dplyr':
## 
##     where
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## Attaching package: 'lme4'
## The following object is masked from 'package:nlme':
## 
##     lmList
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## 3 files sourced from folder "/Users/michaelsieler/Dropbox/Mac (2)/Documents/Sharpton_Lab/Projects_Repository/Rules_of_Life/major-experiment-2023/Code/R/Functions/StartFunctions"
sourceFolder(here::here("Code", "R", "Functions", "HelperFunctions"), T)
## 9 files sourced from folder "/Users/michaelsieler/Dropbox/Mac (2)/Documents/Sharpton_Lab/Projects_Repository/Rules_of_Life/major-experiment-2023/Code/R/Functions/HelperFunctions"
# sourceFolder(here::here("Code", "R", "Functions", "AnalysisScripts"), T)

# Source the calculate_dynamic_metrics.R script to load required functions
source(here::here("Code", "R", "Functions", "AnalysisScripts", "calculate_dynamic_metrics.R"))

Import Data

ps.tmp <- readRDS("/Users/michaelsieler/Dropbox/Mac (2)/Documents/Sharpton_Lab/Projects_Repository/Rules_of_Life/major-experiment-2023/Data/Robjects/pseq_uncleaned_05052025.rds")


ps.cleaned <-
    ps.tmp %>%
        ## Update Metadata
        ps_rename(Time = Timepoint) %>%
        microViz::ps_mutate(
            Treatment = case_when(
                Antibiotics == 0 & Temperature == 0 & Pathogen == 0 ~ "A- T- P-",
                Antibiotics == 0 & Temperature == 0 & Pathogen == 1 ~ "A- T- P+",
                Antibiotics == 1 & Temperature == 0 & Pathogen == 0 ~ "A+ T- P-",
                Antibiotics == 1 & Temperature == 0 & Pathogen == 1 ~ "A+ T- P+",
                Antibiotics == 0 & Temperature == 1 & Pathogen == 0 ~ "A- T+ P-",
                Antibiotics == 0 & Temperature == 1 & Pathogen == 1 ~ "A- T+ P+",
                Antibiotics == 1 & Temperature == 1 & Pathogen == 0 ~ "A+ T+ P-",
                Antibiotics == 1 & Temperature == 1 & Pathogen == 1 ~ "A+ T+ P+",
                TRUE ~ "Unknown"
            ), .after = "Pathogen"
        ) %>%
        microViz::ps_mutate(Sample = fecal.sample.number, .before = 1) %>%
        microViz::ps_mutate(Sample = gsub("^f", "", Sample)) %>%
        microViz::ps_filter(Treatment != "Unknown") %>%
        microViz::ps_mutate(
            History = case_when(
                Antibiotics + Temperature == 0 ~ 0,
                Antibiotics + Temperature == 1 ~ 1,
                Antibiotics + Temperature == 2 ~ 2,
            ), .after = "Treatment"
        ) %>%
        
        ## Additional metadata updates, factorizing metadata
        microViz::ps_mutate(
        # Create treatment code
            treatment_code = case_when(
              Antibiotics == 0 & Temperature == 0 & Pathogen == 0 ~ "Aneg_Tneg_Pneg",
              Antibiotics == 0 & Temperature == 0 & Pathogen == 1 ~ "Aneg_Tneg_Ppos",
              Antibiotics == 1 & Temperature == 0 & Pathogen == 0 ~ "Apos_Tneg_Pneg",
              Antibiotics == 1 & Temperature == 0 & Pathogen == 1 ~ "Apos_Tneg_Ppos",
              Antibiotics == 0 & Temperature == 1 & Pathogen == 0 ~ "Aneg_Tpos_Pneg",
              Antibiotics == 0 & Temperature == 1 & Pathogen == 1 ~ "Aneg_Tpos_Ppos",
              Antibiotics == 1 & Temperature == 1 & Pathogen == 0 ~ "Apos_Tpos_Pneg",
              Antibiotics == 1 & Temperature == 1 & Pathogen == 1 ~ "Apos_Tpos_Ppos"
            ),
            # Create treatment group factor
            treatment_group = case_when(
              Antibiotics == 0 & Temperature == 0 & Pathogen == 1 ~ "Parasite",
              Antibiotics == 1 & Temperature == 0 & Pathogen == 0 ~ "Antibiotics",
              Antibiotics == 1 & Temperature == 0 & Pathogen == 1 ~ "Antibiotics_Parasite",
              Antibiotics == 0 & Temperature == 1 & Pathogen == 0 ~ "Temperature",
              Antibiotics == 0 & Temperature == 1 & Pathogen == 1 ~ "Temperature_Parasite",
              Antibiotics == 1 & Temperature == 1 & Pathogen == 0 ~ "Antibiotics_Temperature",
              Antibiotics == 1 & Temperature == 1 & Pathogen == 1 ~ "Antibiotics_Temperature_Parasite",
              TRUE ~ "Control"
            ),
            # Convert to factor with appropriate levels
            treatment_group = factor(treatment_group, 
                                   levels = c("Control", "Parasite", 
                                              "Antibiotics", "Antibiotics_Parasite",
                                              "Temperature", "Temperature_Parasite",
                                            "Antibiotics_Temperature", "Antibiotics_Temperature_Parasite")
                                   ),
            treatment_code = factor(treatment_code, levels = treatment_order),
            # Create time point factor
            time_point = factor(Time, levels = c(0, 14, 18, 25, 29, 60)),
            # Create pathogen status factor
            pathogen_status = factor(ifelse(Pathogen == 1, "Exposed", "Unexposed"),
                                   levels = c("Unexposed", "Exposed")),
            # Create sex factor
            sex = factor(Sex, levels = c("M", "F"))
            )  %>%
    microViz::ps_mutate(Treatment = factor(Treatment, levels = treatment_order)) %>%
      microViz::ps_mutate(Exp_Type = case_when(
          Treatment %in% c("A- T- P-", "A- T- P+")  ~ "No prior stressor(s)",
          Treatment %in% c("A+ T- P-", "A+ T- P+")  ~ "Antibiotics",
          Treatment %in% c("A- T+ P-", "A- T+ P+") ~ "Temperature",
          Treatment %in% c("A+ T+ P-", "A+ T+ P+") ~ "Combined",
      )) %>%
      microViz::ps_mutate(Exp_Type = factor(Exp_Type, levels = c("No prior stressor(s)", "Antibiotics", "Temperature", "Combined"))) %>%
  # Fix names for taxonomic ranks not identified
  microViz::tax_fix(suffix_rank = "current", anon_unique = T, unknown = NA) %>% 
  # Filter for any samples that contain more than 5000 reads
  microViz::ps_filter(sample_sums(.) > 5000) %>%
  # Any taxa not found in at least 3 samples are removed
  microViz::tax_filter(min_prevalence = 3, undetected = 0) %>%
  # Remove any unwanted reads
  microViz::tax_select(c("Mitochondria", "Chloroplast", "Eukaryota"), deselect = TRUE) %>%
  microViz::tax_select(c("Bacteria, Phylum"), deselect = TRUE) 


# ps.cleaned %>% microViz::samdat_tbl()

Variables

## Diversity Methods

diversity.method <- list()
diversity.method[["alpha"]] <- c("shannon", 
                                 "inverse_simpson", 
                                 "observed")
diversity.method[["beta"]] <- c("bray", "canberra")

## Stats/Plotting Variables

alpha.stats <- list() # save all stat results here
alpha.plots <- list() # save plots here

beta.stats <- list() # save all stat results here
beta.plots <- list() # save plots here

diffAbnd.stats <- list() # save all stat results here
diffAbnd.plots <- list() # save plots here

Calculate Scores

Alpha

## Alpha
ps.cleaned_alpha <-
    ps.cleaned %>% # phyloseq object we'll be using
  microViz::ps_calc_diversity( # calculates various alpha diversity metrics
    rank = "Genus", # What taxonomic rank do you want to calculate diversity metrics at
    index = "shannon", # What diversity metric to use
    varname = "Shannon__Genus", # Column name in sample data
    exp = F # exponentiate the result or not
  ) %>%
  microViz::ps_calc_diversity(
    rank = "Genus",
    index = "gini_simpson",
    varname = "Simpson__Genus" # Column name in sample data
  ) %>%
  microViz::ps_calc_richness( # related to ps_calc_diversity, this calculates richness or "observed" values
    rank = "Genus",
    varname = "Richness__Genus"
  ) %>%
  # ps_calc_diversity.phy(
  #   varname = "Phylogenetic__Genus"
  # ) %>% # Note: This is a helper function I made. You need to adjust function manually to change taxon rank (See: microViz_Helper.R)
  # # ps_calc_diversity.phy() helper function does not properly name the column. Modify the following function as needed 
  # ps_rename("Phylogenetic__Genus" = "phylogenetic_Genus") %>% # Helper function (See: microViz_Helper.R)
  ps_mutate(across(contains("__Genus"), norm_scores, .names = "{.col}_norm")) 
## 
##     lambda     W Shapiro.p.value    A Anderson.p.value
## 420  0.475 0.987       3.787e-05 2.49          2.7e-06
## 
## if (lambda >  0){TRANS = x ^ lambda} 
## if (lambda == 0){TRANS = log(x)} 
## if (lambda <  0){TRANS = -1 * x ^ lambda} 
## 
## 
##     lambda      W Shapiro.p.value     A Anderson.p.value
## 482  2.025 0.9761       2.845e-08 2.862         3.33e-07
## 
## if (lambda >  0){TRANS = x ^ lambda} 
## if (lambda == 0){TRANS = log(x)} 
## if (lambda <  0){TRANS = -1 * x ^ lambda} 
## 
## 
##     lambda      W Shapiro.p.value    A Anderson.p.value
## 393   -0.2 0.9883       0.0001119 1.69        0.0002453
## 
## if (lambda >  0){TRANS = x ^ lambda} 
## if (lambda == 0){TRANS = log(x)} 
## if (lambda <  0){TRANS = -1 * x ^ lambda}
ps.cleaned_alpha %>%
    microViz::samdat_tbl()
## # A tibble: 594 × 37
##    .sample_name Sample Tank.ID  Time Antibiotics Temperature Pathogen Treatment
##    <chr>        <chr>    <dbl> <dbl>       <dbl>       <dbl>    <dbl> <fct>    
##  1 f1           1            1     0           0           0        0 A- T- P- 
##  2 f2           2            1     0           0           0        0 A- T- P- 
##  3 f3           3            1     0           0           0        0 A- T- P- 
##  4 f31          31           2     0           0           0        0 A- T- P- 
##  5 f32          32           2     0           0           0        0 A- T- P- 
##  6 f33          33           2     0           0           0        0 A- T- P- 
##  7 f61          61           3     0           0           0        0 A- T- P- 
##  8 f62          62           3     0           0           0        0 A- T- P- 
##  9 f63          63           3     0           0           0        0 A- T- P- 
## 10 f91          91           4     0           0           0        1 A- T- P+ 
## # ℹ 584 more rows
## # ℹ 29 more variables: History <dbl>, Sex <chr>, Weight.g <dbl>,
## #   Weight.mg <dbl>, Length.mm <dbl>, Body.Size <chr>, Mortality <dbl>,
## #   Total.Worm.Count <dbl>, Flagged <dbl>, notes <chr>,
## #   fecal.sample.number <chr>, gut.sample.number <chr>, adultworm.count <dbl>,
## #   larvalworm.count <dbl>, notes.dissection <chr>, PLATE <chr>, WELL <chr>,
## #   treatment_code <fct>, treatment_group <fct>, time_point <fct>, …
test.data.alpha <-
    ps.cleaned_alpha %>% 
    microViz::samdat_tbl() %>%
    dplyr::relocate(c(".sample_name"), .before = "Antibiotics") %>%
    tidyr::pivot_longer(cols = contains("__Genus"), names_to = "Alpha.Metric", values_to = "Alpha.Score") %>%
    dplyr::mutate(Data.Type = case_when(
        Alpha.Metric %in% c("Shannon__Genus_norm", "Simpson__Genus_norm", "Richness__Genus_norm") ~ "Normalized",
        TRUE ~ "Raw"
    )) %>%
    dplyr::relocate(c("Alpha.Metric", "Alpha.Score", "Data.Type"), .after = "Tank.ID") %>%
    dplyr::relocate(c("Sample", "Treatment", "Time", "Tank.ID"), .before = 1)
    # tidyr::pivot_longer(cols = contains("_norm"), names_to = "Alpha.Metric", values_to = "Alpha.Score") %>%
    # dplyr::filter(Alpha.Metric %in% c("Simpson__Genus_norm", "Simpson__Genus_norm")) %>%
    # tidyr::pivot_longer(cols = contains("__Genus"), names_to = "Raw.Metric", values_to = "Raw.Score") 

Beta

Dist Matrices
#### Beta diversity matrices

##### All

beta.dist.mat <- list()

beta.dist.mat[["All"]] <- # saves the results of this loop to a list
  furrr::future_map(diversity.method[["beta"]], function(x){ # Future_map runs the function in parallel
    ps.cleaned_alpha %>% # is the phyloseq object we are looping over
      microViz::tax_transform(trans = "identity", # tax transform transforms taxa counts
                    rank = ifelse(x == "gunifrac" || x == "wunifrac" || x == "unifrac", "unique", "Genus")) %>% # phylogenetic methods need "unique" ASV counts
      microViz::dist_calc(x, gunifrac_alpha = 0.5) # calculates the distance between samples. gunifrac_alpha = 0.5 weights abundance into its calculation
  }) %>%
  setNames(diversity.method[["beta"]]) # assigns names of the beta methods to the list

beta.dist.mat[["All_60DPE"]] <- # saves the results of this loop to a list
  furrr::future_map(diversity.method[["beta"]], function(x){ # Future_map runs the function in parallel
    ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Time == 60) %>%
      microViz::tax_transform(trans = "identity", # tax transform transforms taxa counts
                    rank = ifelse(x == "gunifrac" || x == "wunifrac" || x == "unifrac", "unique", "Genus")) %>% # phylogenetic methods need "unique" ASV counts
      microViz::dist_calc(x, gunifrac_alpha = 0.5) # calculates the distance between samples. gunifrac_alpha = 0.5 weights abundance into its calculation
  }) %>%
  setNames(diversity.method[["beta"]]) # assigns names of the beta methods to the list

##### Parasite exposure only

beta.dist.mat[["ParasiteExp"]] <- # saves the results of this loop to a list
  furrr::future_map(diversity.method[["beta"]], function(x){ # Future_map runs the function in parallel
    ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Treatment  %in% c(
          "A- T- P-",
          "A- T- P+"#,
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          )) %>%
      microViz::tax_transform(trans = "identity", # tax transform transforms taxa counts
                    rank = ifelse(x == "gunifrac" || x == "wunifrac" || x == "unifrac", "unique", "Genus")) %>% # phylogenetic methods need "unique" ASV counts
      microViz::dist_calc(x, gunifrac_alpha = 0.5) # calculates the distance between samples. gunifrac_alpha = 0.5 weights abundance into its calculation
  }) %>%
  setNames(diversity.method[["beta"]]) # assigns names of the beta methods to the list

beta.dist.mat[["ParasiteExp_60DPE"]] <- # saves the results of this loop to a list
  furrr::future_map(diversity.method[["beta"]], function(x){ # Future_map runs the function in parallel
    ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Treatment  %in% c(
          "A- T- P-",
          "A- T- P+"#,
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          )) %>%
          microViz::ps_filter(Time == 60) %>%
      microViz::tax_transform(trans = "identity", # tax transform transforms taxa counts
                    rank = ifelse(x == "gunifrac" || x == "wunifrac" || x == "unifrac", "unique", "Genus")) %>% # phylogenetic methods need "unique" ASV counts
      microViz::dist_calc(x, gunifrac_alpha = 0.5) # calculates the distance between samples. gunifrac_alpha = 0.5 weights abundance into its calculation
  }) %>%
  setNames(diversity.method[["beta"]]) # assigns names of the beta methods to the list

##### Prior stressors and parasite exposure

beta.dist.mat[["PriorStressParaExp"]] <- # saves the results of this loop to a list
  furrr::future_map(diversity.method[["beta"]], function(x){ # Future_map runs the function in parallel
    ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Treatment  %in% c(
          # "A- T- P-",
          "A- T- P+",
          # "A+ T- P-",
          "A+ T- P+",
          # "A- T+ P-",
          "A- T+ P+",
          # "A+ T+ P-",
          "A+ T+ P+"
          )) %>%
      microViz::tax_transform(trans = "identity", # tax transform transforms taxa counts
                    rank = ifelse(x == "gunifrac" || x == "wunifrac" || x == "unifrac", "unique", "Genus")) %>% # phylogenetic methods need "unique" ASV counts
      microViz::dist_calc(x, gunifrac_alpha = 0.5) # calculates the distance between samples. gunifrac_alpha = 0.5 weights abundance into its calculation
  }) %>%
  setNames(diversity.method[["beta"]]) # assigns names of the beta methods to the list

beta.dist.mat[["PriorStressParaExp_60DPE"]] <- # saves the results of this loop to a list
  furrr::future_map(diversity.method[["beta"]], function(x){ # Future_map runs the function in parallel
    ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Treatment  %in% c(
          # "A- T- P-",
          "A- T- P+",
          # "A+ T- P-",
          "A+ T- P+",
          # "A- T+ P-",
          "A- T+ P+",
          # "A+ T+ P-",
          "A+ T+ P+"
          )) %>%
          microViz::ps_filter(Time == 60) %>%
      microViz::tax_transform(trans = "identity", # tax transform transforms taxa counts
                    rank = ifelse(x == "gunifrac" || x == "wunifrac" || x == "unifrac", "unique", "Genus")) %>% # phylogenetic methods need "unique" ASV counts
      microViz::dist_calc(x, gunifrac_alpha = 0.5) # calculates the distance between samples. gunifrac_alpha = 0.5 weights abundance into its calculation
  }) %>%
  setNames(diversity.method[["beta"]]) # assigns names of the beta methods to the list

##### Prior stressors and NO parasite exposure

beta.dist.mat[["PriorStressNoParaExp"]] <- # saves the results of this loop to a list
  furrr::future_map(diversity.method[["beta"]], function(x){ # Future_map runs the function in parallel
    ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Treatment  %in% c(
          "A- T- P-",
          # "A- T- P+",
          "A+ T- P-",
          # "A+ T- P+",
          "A- T+ P-",
          # "A- T+ P+",
          "A+ T+ P-"#,
          # "A+ T+ P+"
          )) %>%
      microViz::tax_transform(trans = "identity", # tax transform transforms taxa counts
                    rank = ifelse(x == "gunifrac" || x == "wunifrac" || x == "unifrac", "unique", "Genus")) %>% # phylogenetic methods need "unique" ASV counts
      microViz::dist_calc(x, gunifrac_alpha = 0.5) # calculates the distance between samples. gunifrac_alpha = 0.5 weights abundance into its calculation
  }) %>%
  setNames(diversity.method[["beta"]]) # assigns names of the beta methods to the list

beta.dist.mat[["PriorStressNoParaExp_60DPE"]] <- # saves the results of this loop to a list
  furrr::future_map(diversity.method[["beta"]], function(x){ # Future_map runs the function in parallel
    ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Treatment  %in% c(
          "A- T- P-",
          # "A- T- P+",
          "A+ T- P-",
          # "A+ T- P+",
          "A- T+ P-",
          # "A- T+ P+",
          "A+ T+ P-"#,
          # "A+ T+ P+"
          )) %>%
          microViz::ps_filter(Time == 60) %>%
      microViz::tax_transform(trans = "identity", # tax transform transforms taxa counts
                    rank = ifelse(x == "gunifrac" || x == "wunifrac" || x == "unifrac", "unique", "Genus")) %>% # phylogenetic methods need "unique" ASV counts
      microViz::dist_calc(x, gunifrac_alpha = 0.5) # calculates the distance between samples. gunifrac_alpha = 0.5 weights abundance into its calculation
  }) %>%
  setNames(diversity.method[["beta"]]) # assigns names of the beta methods to the list
Personalization
#### Beta Personalization Scores

data.beta.list <- list()

# Save dataframe
data.beta.list[["All"]] <-
  
  ps.cleaned_alpha %>% 
  
  # Calculate Distance matrix
  microViz::dist_calc("bray") %>% 
  
  # Get distance matric
  microViz::dist_get() %>%
  # Pat Schloss methods to view beta diversity matrix data
  # Source: https://youtu.be/Gm-kg3TuML4?si=nYLLAxm7uv8aYG27&t=250
  
  # Convert into a matrix object
  as.matrix() %>%
  
  # Convert matrix into a tibble
  tibble::as_tibble(rownames = "Sample_a") %>%
  
  # Convert dataframe into "long" format
  tidyr::pivot_longer(-Sample_a, names_to = "Sample_b", values_to = "Dist") %>%
  
  # Remove any rows where names equal eachother
  dplyr::filter(Sample_a != Sample_b) %>%
  
  # Merge sample dataframe
  dplyr::right_join((ps.cleaned_alpha %>% 
                       microViz::samdat_tbl() %>% ungroup() %>%  arrange(Time, Treatment, Tank.ID) %>% 
                       dplyr::select(!contains("_Genus")) %>%
                       dplyr::select(fecal.sample.number, Time, Treatment, Tank.ID)
  ), 
  by = c("Sample_a" = "fecal.sample.number")) %>% 
  dplyr::right_join((ps.cleaned_alpha %>% 
                       microViz::samdat_tbl() %>% 
                       dplyr::ungroup() %>% 
                       dplyr::arrange(Time, Treatment, Tank.ID) %>% 
                       dplyr::select(!contains("_Genus")) %>%
                       dplyr::select(fecal.sample.number, Time, Treatment, Tank.ID)
  ), 
  by = c("Sample_b" = "fecal.sample.number"), suffix = c("_a", "_b")) %>% # add suffix to differentiate between sample metadata
  # Move columns around
  dplyr::select(order(colnames(.))) %>%
  dplyr::relocate(dplyr::any_of(c("Sample_a", "Sample_b", "Dist")), .before = 1) %>%
  dplyr::relocate(dplyr::contains(c("Tank")), .after = dplyr::last_col()) %>%
  dplyr::ungroup() %>%
  
  # Calculate the median distance of samples to all other samples (Measure of "personalization")
  dplyr::group_by(Sample_a, Treatment_a, Time_a, Tank.ID_a) %>%
  dplyr::summarize(Beta.Score = median(Dist), #*100,  
                   Beta.Score.Mean = mean(Dist), #*100, 
                   count = dplyr::n()) %>% 
  dplyr::ungroup() %>%
  
  # Rename Columns
  dplyr::rename(fecal.sample.number = Sample_a,
                Treatment = Treatment_a,
                
                Time = Time_a,
                Tank.ID = Tank.ID_a) %>%
  dplyr::right_join(ps.cleaned_alpha %>%
                      microViz::samdat_tbl() %>% 
                      dplyr::select(!contains("_Genus")),
                    by = c("fecal.sample.number", "Time", "Treatment","Tank.ID")) %>%
  dplyr::ungroup() %>%
  
  # Normalize diversity scores
  dplyr::mutate(dplyr::across(dplyr::contains("Beta.Score"), norm_scores, .names = "{.col}_norm"), .after = "Beta.Score") %>%
  
  # Create a column to save what type of beta metric was calculated
  dplyr::mutate(Beta.Metric = "Bray-Curtis", .before = "Beta.Score") %>%
  dplyr::mutate(Sample = gsub("^f", "", fecal.sample.number)) %>%
    dplyr::relocate(Sample, .before = 1) %>%
    dplyr::relocate(fecal.sample.number, .after = 36)
## `summarise()` has grouped output by 'Sample_a', 'Treatment_a', 'Time_a'. You
## can override using the `.groups` argument.
## 
##     lambda      W Shapiro.p.value    A Anderson.p.value
## 291  -2.75 0.9882       0.0001004 1.05         0.009232
## 
## if (lambda >  0){TRANS = x ^ lambda} 
## if (lambda == 0){TRANS = log(x)} 
## if (lambda <  0){TRANS = -1 * x ^ lambda} 
## 
## 
##     lambda      W Shapiro.p.value     A Anderson.p.value
## 235  -4.15 0.9883       0.0001067 1.004          0.01196
## 
## if (lambda >  0){TRANS = x ^ lambda} 
## if (lambda == 0){TRANS = log(x)} 
## if (lambda <  0){TRANS = -1 * x ^ lambda}
data.beta.list[["All"]] <- data.beta.list[["All"]] %>%
dplyr::full_join(., 
                   ps.cleaned_alpha %>% 
                     # tax_agg("Genus") %>%
                     microViz::dist_calc("canberra") %>% 
                     
                     # Get distance matric
    microViz::dist_get() %>%
    # Pat Schloss methods to view beta diversity matrix data
    # Source: https://youtu.be/Gm-kg3TuML4?si=nYLLAxm7uv8aYG27&t=250
    
    # Convert into a matrix object
    as.matrix() %>%
    
    # Convert matrix into a tibble
    tibble::as_tibble(rownames = "Sample_a") %>%
    
    # Convert dataframe into "long" format
    tidyr::pivot_longer(-Sample_a, names_to = "Sample_b", values_to = "Dist") %>%
    
    # Remove any rows where names equal eachother
    dplyr::filter(Sample_a != Sample_b) %>%
    
    # Merge sample dataframe
    dplyr::right_join((ps.cleaned_alpha %>% 
                       microViz::samdat_tbl() %>% ungroup() %>%  arrange(Time, Treatment, Tank.ID) %>% 
                       dplyr::select(!contains("_Genus")) %>%
                       dplyr::select(fecal.sample.number, Time, Treatment, Tank.ID)
    ), 
    by = c("Sample_a" = "fecal.sample.number")) %>% 
    dplyr::right_join((ps.cleaned_alpha %>% 
                       microViz::samdat_tbl() %>% 
                       dplyr::ungroup() %>% 
                       dplyr::arrange(Time, Treatment, Tank.ID) %>% 
                       dplyr::select(!contains("_Genus")) %>%
                       dplyr::select(fecal.sample.number, Time, Treatment, Tank.ID)
    ), 
    by = c("Sample_b" = "fecal.sample.number"), suffix = c("_a", "_b")) %>% # add suffix to differentiate between sample metadata
    # Move columns around
    dplyr::select(order(colnames(.))) %>%
    dplyr::relocate(dplyr::any_of(c("Sample_a", "Sample_b", "Dist")), .before = 1) %>%
    dplyr::relocate(dplyr::contains(c("Tank")), .after = dplyr::last_col()) %>%
    dplyr::ungroup() %>%
    
    # Calculate the median distance of samples to all other samples (Measure of "personalization")
    dplyr::group_by(Sample_a, Treatment_a, Time_a, Tank.ID_a) %>%
    dplyr::summarize(Beta.Score = median(Dist), #*100,  
                   Beta.Score.Mean = mean(Dist), #*100, 
                   count = dplyr::n()) %>% 
    dplyr::ungroup() %>%
    
    # Rename Columns
    dplyr::rename(fecal.sample.number = Sample_a,
                Treatment = Treatment_a,
                
                Time = Time_a,
                Tank.ID = Tank.ID_a) %>%
    dplyr::right_join(ps.cleaned_alpha %>%
                      microViz::samdat_tbl() %>% 
                      dplyr::select(!contains("_Genus")),
                    by = c("fecal.sample.number", "Time", "Treatment","Tank.ID")) %>%
    dplyr::ungroup() %>%
                     
    # Normalize diversity scores
    dplyr::mutate(dplyr::across(contains("Beta.Score"), norm_scores, .names = "{.col}_norm"), .after = "Beta.Score") %>%
    
    # Create a column to save what type of beta metric was calculated
    dplyr::mutate(Beta.Metric = "Canberra", .before = "Beta.Score")
  )
## `summarise()` has grouped output by 'Sample_a', 'Treatment_a', 'Time_a'. You
## can override using the `.groups` argument.
## 
##    lambda      W Shapiro.p.value      A Anderson.p.value
## 76 -8.125 0.9925        0.004258 0.6477          0.09061
## 
## if (lambda >  0){TRANS = x ^ lambda} 
## if (lambda == 0){TRANS = log(x)} 
## if (lambda <  0){TRANS = -1 * x ^ lambda} 
## 
## 
##    lambda      W Shapiro.p.value      A Anderson.p.value
## 51  -8.75 0.9948         0.04029 0.4685           0.2483
## 
## if (lambda >  0){TRANS = x ^ lambda} 
## if (lambda == 0){TRANS = log(x)} 
## if (lambda <  0){TRANS = -1 * x ^ lambda}
## Joining with `by = join_by(Sample, Treatment, Time, Tank.ID, Beta.Metric,
## Beta.Score, Beta.Score_norm, Beta.Score.Mean_norm, Beta.Score.Mean, count,
## .sample_name, Antibiotics, Temperature, Pathogen, History, Sex, Weight.g,
## Weight.mg, Length.mm, Body.Size, Mortality, Total.Worm.Count, Flagged, notes,
## gut.sample.number, adultworm.count, larvalworm.count, notes.dissection, PLATE,
## WELL, treatment_code, treatment_group, time_point, pathogen_status, sex,
## fecal.sample.number, Exp_Type)`
# Prepare beta diversity data
test.data.beta <-
    data.beta.list[["All"]] %>%
        tidyr::pivot_longer(cols = c("Beta.Score", "Beta.Score_norm"), names_to = "Beta.Metric.Type", values_to = "Beta.Score") %>%
        dplyr::relocate(c("Beta.Metric.Type", "Beta.Score"), .after = "Beta.Metric") %>%
        dplyr::mutate(Data.Type = case_when(
            Beta.Metric.Type %in% c("   Beta.Score", "Beta.Score_norm") ~ "Normalized",
            TRUE ~ "Raw"
        ), .after = "Beta.Score") %>%
        # Remove the temporary Beta.Metric.Type column
        dplyr::select(-c(Beta.Metric.Type, count)) %>%
        dplyr::mutate(Beta.Metric = ifelse(Data.Type == "Normalized", paste0(Beta.Metric, "_norm"), Beta.Metric)) %>%
        dplyr::arrange(Sample, Time, Tank.ID, Data.Type)

Calculate SADMMs

# Calculate all metrics at once using the combined function
test.alpha.all_metrics <- calculate_all_dynamic_metrics(
  data = test.data.alpha,
  metric_col = "Alpha.Metric",
  score_col = "Alpha.Score",
  ref_group = NULL,  # Using all groups as reference
  time_col = "Time",
  group_col = "Treatment",
  use_mean = TRUE
)
## Joining with `by = join_by(Sample, Time, Treatment, Tank.ID, Alpha.Metric,
## Data.Type, Displacement)`
## Joining with `by = join_by(Sample, Time, Treatment, Tank.ID, Alpha.Metric,
## Data.Type)`
# Calculate all metrics at once using the combined function
test.beta.all_metrics <- calculate_all_dynamic_metrics(
  data = test.data.beta,
  metric_col = "Beta.Metric",
  score_col = "Beta.Score",
  ref_group = NULL,
  time_col = "Time",
  group_col = "Treatment",
  use_mean = TRUE
)
## Joining with `by = join_by(Sample, Time, Treatment, Tank.ID, Beta.Metric,
## Data.Type, Displacement)`
## Joining with `by = join_by(Sample, Time, Treatment, Tank.ID, Beta.Metric,
## Data.Type)`
# Merge Displacement data with Alpha.Metric and Alpha.Score
data.alpha.merged <-
    test.alpha.all_metrics$volatility %>%
        # dplyr::filter()
        # First pivot the SADMM metrics
        # dplyr::group_by(Data.Type) %>%
        tidyr::pivot_longer(
            cols = c("Alpha.Score", "Ref.Score", "Displacement", "Velocity", "Acceleration", "Volatility"),
            names_to = "METRIC",
            values_to = "SCORE"
        ) %>%
        dplyr::mutate(METRIC = factor(METRIC, levels = c("Alpha.Score", "Ref.Score", "Displacement", "Velocity", "Acceleration", "Volatility")),
                      Data.Type = factor(Data.Type, levels = c("Raw", "Normalized")))



# Merge Displacement data with Alpha.Metric and Alpha.Score
data.beta.merged <-
    test.beta.all_metrics$volatility %>%
        # dplyr::filter()
        # First pivot the SADMM metrics
        # dplyr::group_by(Data.Type) %>%
        tidyr::pivot_longer(
            cols = c("Beta.Score", "Ref.Score", "Displacement", "Velocity", "Acceleration", "Volatility"),
            names_to = "METRIC",
            values_to = "SCORE"
        ) %>%
        dplyr::mutate(METRIC = factor(METRIC, levels = c("Beta.Score", "Ref.Score", "Displacement", "Velocity", "Acceleration", "Volatility")),
                      Data.Type = factor(Data.Type, levels = c("Raw", "Normalized")))

00) All

Alpha

All

Displacement
Plot

Code
# Set the subsection for this analysis
tmp.resSubSection <- "All"

# GLM Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]] <- 
  test.alpha.all_metrics$volatility %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models(formula_str = "Displacement ~ Treatment", 
                   alpha_metric_col = "Alpha.Metric", 
                   alpha_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["ANOVA"]] <-
  run_glm_anova(alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]])

# Tukey Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["Tukey"]] <-
  test.alpha.all_metrics$volatility %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm(., "Displacement", "Alpha.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
Tables
GLM Results
glm(Displacement ~ Treatment); All treatments
term estimate std.error statistic p.value p.adj.sig
Shannon
(Intercept) 0.078 0.019 3.995 <0.001 ****
A- T- P+ −0.045 0.027 −1.627 0.104 ns
A+ T- P- −0.051 0.028 −1.815 0.070 ns
A+ T- P+ −0.071 0.028 −2.584 0.010 **
A- T+ P- −0.004 0.028 −0.151 ≥0.25 ns
A- T+ P+ −0.092 0.028 −3.302 0.001 **
A+ T+ P- −0.113 0.030 −3.799 <0.001 ***
A+ T+ P+ −0.061 0.029 −2.149 0.032 *
Simpson
(Intercept) 0.042 0.028 1.527 0.127 ns
A- T- P+ −0.049 0.039 −1.238 0.216 ns
A+ T- P- −0.034 0.040 −0.841 ≥0.25 ns
A+ T- P+ −0.085 0.039 −2.165 0.031 *
A- T+ P- 0.035 0.041 0.874 ≥0.25 ns
A- T+ P+ −0.154 0.040 −3.867 <0.001 ***
A+ T+ P- −0.181 0.042 −4.270 <0.001 ****
A+ T+ P+ −0.091 0.041 −2.232 0.026 *
Richness
(Intercept) 0.221 0.024 9.249 <0.001 ****
A- T- P+ 0.009 0.034 0.265 ≥0.25 ns
A+ T- P- −0.049 0.035 −1.430 0.153 ns
A+ T- P+ −0.045 0.034 −1.325 0.186 ns
A- T+ P- −0.016 0.035 −0.467 ≥0.25 ns
A- T+ P+ 0.016 0.034 0.456 ≥0.25 ns
A+ T+ P- 0.029 0.037 0.783 ≥0.25 ns
A+ T+ P+ −0.054 0.035 −1.533 0.126 ns
ANOVA of GLM
ANOVA(GLM(Displacement ~ Treatment), type = 2); All treatments
term statistic df p.value sig
Shannon
Treatment 25.086 7.000 <0.001 ***
Simpson
Treatment 42.565 7.000 <0.001 ****
Richness
Treatment 11.161 7.000 0.132 ns
Pairwise Tukey's HSD, p.adj: Dunnett
Tukey(Displacement ~ Treatment); All treatments
term .y. group1 group2 estimate std.error statistic adj.p.value Variable
Shannon
Treatment Alpha.Score A- T- P+ A- T- P- −0.045 0.027 −1.627 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P- A- T- P- −0.051 0.028 −1.815 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P+ A- T- P- −0.071 0.028 −2.584 0.161 Treatment
Treatment Alpha.Score A- T+ P- A- T- P- −0.004 0.028 −0.151 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A- T- P- −0.092 0.028 −3.302 0.022 Treatment
Treatment Alpha.Score A+ T+ P- A- T- P- −0.113 0.030 −3.799 0.004 Treatment
Treatment Alpha.Score A+ T+ P+ A- T- P- −0.061 0.029 −2.149 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P- A- T- P+ −0.006 0.028 −0.220 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P+ A- T- P+ −0.027 0.028 −0.968 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A- T- P+ 0.040 0.028 1.421 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A- T- P+ −0.047 0.028 −1.695 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A- T- P+ −0.068 0.030 −2.298 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A- T- P+ −0.017 0.029 −0.582 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P+ A+ T- P- −0.021 0.028 −0.730 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A+ T- P- 0.047 0.029 1.608 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A+ T- P- −0.041 0.028 −1.446 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A+ T- P- −0.062 0.030 −2.056 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A+ T- P- −0.010 0.029 −0.360 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A+ T- P+ 0.067 0.029 2.348 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A+ T- P+ −0.020 0.028 −0.729 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A+ T- P+ −0.042 0.030 −1.392 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A+ T- P+ 0.010 0.029 0.353 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A- T+ P- −0.088 0.029 −3.044 0.048 Treatment
Treatment Alpha.Score A+ T+ P- A- T+ P- −0.109 0.031 −3.549 0.009 Treatment
Treatment Alpha.Score A+ T+ P+ A- T+ P- −0.057 0.029 −1.936 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A- T+ P+ −0.021 0.030 −0.706 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A- T+ P+ 0.031 0.029 1.058 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A+ T+ P- 0.052 0.031 1.684 ≥0.25 Treatment
Simpson
Treatment Alpha.Score A- T- P+ A- T- P- −0.049 0.039 −1.238 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P- A- T- P- −0.034 0.040 −0.841 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P+ A- T- P- −0.085 0.039 −2.165 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A- T- P- 0.035 0.041 0.874 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A- T- P- −0.154 0.040 −3.867 0.003 Treatment
Treatment Alpha.Score A+ T+ P- A- T- P- −0.181 0.042 −4.270 <0.001 Treatment
Treatment Alpha.Score A+ T+ P+ A- T- P- −0.091 0.041 −2.232 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P- A- T- P+ 0.015 0.040 0.372 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P+ A- T- P+ −0.037 0.039 −0.935 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A- T- P+ 0.084 0.041 2.070 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A- T- P+ −0.105 0.040 −2.645 0.140 Treatment
Treatment Alpha.Score A+ T+ P- A- T- P+ −0.133 0.042 −3.128 0.038 Treatment
Treatment Alpha.Score A+ T+ P+ A- T- P+ −0.042 0.041 −1.041 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P+ A+ T- P- −0.052 0.040 −1.287 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A+ T- P- 0.069 0.041 1.672 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A+ T- P- −0.120 0.040 −2.963 0.060 Treatment
Treatment Alpha.Score A+ T+ P- A+ T- P- −0.148 0.043 −3.421 0.015 Treatment
Treatment Alpha.Score A+ T+ P+ A+ T- P- −0.057 0.041 −1.381 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A+ T- P+ 0.121 0.041 2.962 0.061 Treatment
Treatment Alpha.Score A- T+ P+ A+ T- P+ −0.068 0.040 −1.706 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A+ T- P+ −0.096 0.043 −2.248 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A+ T- P+ −0.005 0.041 −0.134 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A- T+ P- −0.189 0.041 −4.604 <0.001 Treatment
Treatment Alpha.Score A+ T+ P- A- T+ P- −0.217 0.044 −4.957 <0.001 Treatment
Treatment Alpha.Score A+ T+ P+ A- T+ P- −0.126 0.042 −3.006 0.053 Treatment
Treatment Alpha.Score A+ T+ P- A- T+ P+ −0.028 0.043 −0.649 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A- T+ P+ 0.063 0.041 1.521 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A+ T+ P- 0.091 0.044 2.063 ≥0.25 Treatment
Richness
Treatment Alpha.Score A- T- P+ A- T- P- 0.009 0.034 0.265 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P- A- T- P- −0.049 0.035 −1.430 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P+ A- T- P- −0.045 0.034 −1.325 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A- T- P- −0.016 0.035 −0.467 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A- T- P- 0.016 0.034 0.456 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A- T- P- 0.029 0.037 0.783 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A- T- P- −0.054 0.035 −1.533 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P- A- T- P+ −0.058 0.035 −1.691 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P+ A- T- P+ −0.054 0.034 −1.589 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A- T- P+ −0.025 0.035 −0.724 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A- T- P+ 0.007 0.034 0.194 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A- T- P+ 0.020 0.037 0.538 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A- T- P+ −0.063 0.035 −1.789 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P+ A+ T- P- 0.004 0.035 0.122 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A+ T- P- 0.033 0.036 0.925 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A+ T- P- 0.065 0.035 1.860 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A+ T- P- 0.078 0.037 2.093 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A+ T- P- −0.005 0.036 −0.126 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A+ T- P+ 0.029 0.035 0.816 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A+ T- P+ 0.061 0.035 1.761 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A+ T- P+ 0.074 0.037 2.002 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A+ T- P+ −0.009 0.035 −0.248 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A- T+ P- 0.032 0.035 0.903 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A- T+ P- 0.045 0.038 1.193 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A- T+ P- −0.038 0.036 −1.034 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A- T+ P+ 0.013 0.037 0.353 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A- T+ P+ −0.070 0.036 −1.955 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A+ T+ P- −0.083 0.038 −2.180 ≥0.25 Treatment

60 DPE

Displacement
Plot

Code
# Set the subsection for this analysis
tmp.resSubSection <- "All_60DPE"

# GLM Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]] <- 
  test.alpha.all_metrics$displacement %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models(formula_str = "Displacement ~ Treatment", 
                   alpha_metric_col = "Alpha.Metric", 
                   alpha_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["ANOVA"]] <-
  run_glm_anova(alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]])

# Tukey Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["Tukey"]] <-
  test.alpha.all_metrics$displacement %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm(., "Displacement", "Alpha.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
Tables
GLM Results - 60 DPE
glm(Displacement ~ Treatment); All treatments at 60 DPE
term estimate std.error statistic p.value p.adj.sig
Shannon
(Intercept) 0.091 0.037 2.457 0.015 *
A- T- P+ −0.082 0.052 −1.590 0.113 ns
A+ T- P- −0.034 0.055 −0.617 ≥0.25 ns
A+ T- P+ −0.083 0.052 −1.590 0.113 ns
A- T+ P- −0.041 0.057 −0.719 ≥0.25 ns
A- T+ P+ −0.159 0.054 −2.965 0.003 **
A+ T+ P- −0.191 0.068 −2.815 0.005 **
A+ T+ P+ −0.079 0.058 −1.364 0.174 ns
Simpson
(Intercept) 0.008 0.049 0.162 ≥0.25 ns
A- T- P+ −0.070 0.069 −1.011 ≥0.25 ns
A+ T- P- 0.023 0.073 0.311 ≥0.25 ns
A+ T- P+ −0.070 0.070 −1.002 ≥0.25 ns
A- T+ P- 0.030 0.076 0.399 ≥0.25 ns
A- T+ P+ −0.249 0.072 −3.483 <0.001 ***
A+ T+ P- −0.332 0.091 −3.665 <0.001 ***
A+ T+ P+ −0.106 0.077 −1.387 0.167 ns
Richness
(Intercept) 0.338 0.032 10.566 <0.001 ****
A- T- P+ 0.025 0.045 0.552 ≥0.25 ns
A+ T- P- −0.064 0.047 −1.349 0.179 ns
A+ T- P+ −0.108 0.046 −2.363 0.019 *
A- T+ P- −0.089 0.049 −1.804 0.073 ns
A- T+ P+ 0.022 0.047 0.462 ≥0.25 ns
A+ T+ P- 0.067 0.059 1.141 ≥0.25 ns
A+ T+ P+ −0.028 0.050 −0.560 ≥0.25 ns
ANOVA of GLM - 60 DPE
ANOVA(GLM(Displacement ~ Treatment), type = 2); All treatments at 60 DPE
term statistic df p.value sig
Shannon
Treatment 14.694 7.000 0.040 *
Simpson
Treatment 31.743 7.000 <0.001 ****
Richness
Treatment 19.090 7.000 0.008 **
Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); All treatments at 60 DPE
term .y. group1 group2 estimate std.error statistic adj.p.value Variable
Shannon
Treatment Alpha.Score A- T- P+ A- T- P- −0.082 0.052 −1.590 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P- A- T- P- −0.034 0.055 −0.617 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P+ A- T- P- −0.083 0.052 −1.590 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A- T- P- −0.041 0.057 −0.719 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A- T- P- −0.159 0.054 −2.965 0.059 Treatment
Treatment Alpha.Score A+ T+ P- A- T- P- −0.191 0.068 −2.815 0.089 Treatment
Treatment Alpha.Score A+ T+ P+ A- T- P- −0.079 0.058 −1.364 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P- A- T- P+ 0.049 0.054 0.894 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P+ A- T- P+ −0.001 0.052 −0.022 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A- T- P+ 0.041 0.057 0.731 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A- T- P+ −0.077 0.053 −1.442 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A- T- P+ −0.109 0.068 −1.611 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A- T- P+ 0.004 0.057 0.066 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P+ A+ T- P- −0.050 0.055 −0.904 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A+ T- P- −0.007 0.059 −0.122 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A+ T- P- −0.126 0.056 −2.234 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A+ T- P- −0.158 0.070 −2.254 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A+ T- P- −0.045 0.060 −0.748 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A+ T- P+ 0.043 0.057 0.743 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A+ T- P+ −0.076 0.054 −1.402 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A+ T- P+ −0.108 0.068 −1.581 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A+ T- P+ 0.005 0.058 0.085 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A- T+ P- −0.118 0.058 −2.027 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A- T+ P- −0.150 0.072 −2.097 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A- T+ P- −0.038 0.062 −0.607 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A- T+ P+ −0.032 0.069 −0.463 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A- T+ P+ 0.081 0.059 1.368 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A+ T+ P- 0.113 0.072 1.562 ≥0.25 Treatment
Simpson
Treatment Alpha.Score A- T- P+ A- T- P- −0.070 0.069 −1.011 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P- A- T- P- 0.023 0.073 0.311 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P+ A- T- P- −0.070 0.070 −1.002 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A- T- P- 0.030 0.076 0.399 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A- T- P- −0.249 0.072 −3.483 0.012 Treatment
Treatment Alpha.Score A+ T+ P- A- T- P- −0.332 0.091 −3.665 0.006 Treatment
Treatment Alpha.Score A+ T+ P+ A- T- P- −0.106 0.077 −1.387 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P- A- T- P+ 0.092 0.072 1.276 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P+ A- T- P+ 0.000 0.070 −0.005 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A- T- P+ 0.100 0.075 1.326 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A- T- P+ −0.180 0.071 −2.525 0.182 Treatment
Treatment Alpha.Score A+ T+ P- A- T- P+ −0.262 0.090 −2.906 0.070 Treatment
Treatment Alpha.Score A+ T+ P+ A- T- P+ −0.037 0.076 −0.480 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P+ A+ T- P- −0.093 0.073 −1.265 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A+ T- P- 0.008 0.079 0.096 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A+ T- P- −0.272 0.075 −3.632 0.007 Treatment
Treatment Alpha.Score A+ T+ P- A+ T- P- −0.355 0.093 −3.804 0.004 Treatment
Treatment Alpha.Score A+ T+ P+ A+ T- P- −0.129 0.080 −1.617 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A+ T- P+ 0.100 0.076 1.315 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A+ T- P+ −0.179 0.072 −2.488 0.198 Treatment
Treatment Alpha.Score A+ T+ P- A+ T- P+ −0.262 0.091 −2.879 0.075 Treatment
Treatment Alpha.Score A+ T+ P+ A+ T- P+ −0.036 0.077 −0.470 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A- T+ P- −0.280 0.078 −3.594 0.007 Treatment
Treatment Alpha.Score A+ T+ P- A- T+ P- −0.362 0.096 −3.790 0.004 Treatment
Treatment Alpha.Score A+ T+ P+ A- T+ P- −0.137 0.083 −1.655 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A- T+ P+ −0.083 0.092 −0.895 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A- T+ P+ 0.143 0.079 1.818 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A+ T+ P- 0.226 0.096 2.343 ≥0.25 Treatment
Richness
Treatment Alpha.Score A- T- P+ A- T- P- 0.025 0.045 0.552 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P- A- T- P- −0.064 0.047 −1.349 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P+ A- T- P- −0.108 0.046 −2.363 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A- T- P- −0.089 0.049 −1.804 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A- T- P- 0.022 0.047 0.462 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A- T- P- 0.067 0.059 1.141 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A- T- P- −0.028 0.050 −0.560 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P- A- T- P+ −0.089 0.047 −1.884 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P+ A- T- P+ −0.132 0.045 −2.927 0.066 Treatment
Treatment Alpha.Score A- T+ P- A- T- P+ −0.114 0.049 −2.320 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A- T- P+ −0.003 0.046 −0.070 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A- T- P+ 0.043 0.059 0.724 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A- T- P+ −0.053 0.050 −1.062 ≥0.25 Treatment
Treatment Alpha.Score A+ T- P+ A+ T- P- −0.044 0.048 −0.913 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A+ T- P- −0.025 0.051 −0.488 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A+ T- P- 0.086 0.049 1.755 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A+ T- P- 0.131 0.061 2.164 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A+ T- P- 0.036 0.052 0.694 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A+ T- P+ 0.019 0.050 0.373 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A+ T- P+ 0.129 0.047 2.752 0.106 Treatment
Treatment Alpha.Score A+ T+ P- A+ T- P+ 0.175 0.059 2.954 0.061 Treatment
Treatment Alpha.Score A+ T+ P+ A+ T- P+ 0.080 0.050 1.585 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A- T+ P- 0.111 0.051 2.184 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A- T+ P- 0.156 0.062 2.514 0.187 Treatment
Treatment Alpha.Score A+ T+ P+ A- T+ P- 0.061 0.054 1.137 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A- T+ P+ 0.046 0.060 0.762 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A- T+ P+ −0.050 0.051 −0.967 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A+ T+ P- −0.095 0.063 −1.520 ≥0.25 Treatment
Significant Results:Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); All treatments at 60 DPE
term .y. group1 group2 estimate std.error statistic adj.p.value Variable
Simpson
Treatment Alpha.Score A- T+ P+ A- T- P- −0.249 0.072 −3.483 0.012 Treatment
Treatment Alpha.Score A+ T+ P- A- T- P- −0.332 0.091 −3.665 0.006 Treatment
Treatment Alpha.Score A- T+ P+ A+ T- P- −0.272 0.075 −3.632 0.007 Treatment
Treatment Alpha.Score A+ T+ P- A+ T- P- −0.355 0.093 −3.804 0.004 Treatment
Treatment Alpha.Score A- T+ P+ A- T+ P- −0.280 0.078 −3.594 0.007 Treatment
Treatment Alpha.Score A+ T+ P- A- T+ P- −0.362 0.096 −3.790 0.004 Treatment

Beta

All

Displacement
Plot

Code
# Set the subsection for this analysis
tmp.resSubSection <- "All"

# GLM Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]] <- 
  test.beta.all_metrics$volatility %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models_beta(formula_str = "Displacement ~ Treatment", 
                   beta_metric_col = "Beta.Metric", 
                   beta_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["ANOVA"]] <-
  run_glm_anova_beta(beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]])

# Tukey Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["Tukey"]] <-
  test.beta.all_metrics$volatility %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm_beta(., "Displacement", "Beta.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
Tables
GLM Results - Beta
glm(Displacement ~ Treatment); All treatments
term estimate std.error statistic p.value p.adj.sig
Bray-Curtis_norm
(Intercept) −0.061 0.027 −2.266 0.024 *
A- T- P+ 0.142 0.038 3.716 <0.001 ***
A+ T- P- 0.016 0.039 0.412 ≥0.25 ns
A+ T- P+ 0.002 0.039 0.040 ≥0.25 ns
A- T+ P- 0.081 0.040 2.048 0.041 *
A- T+ P+ −0.163 0.039 −4.216 <0.001 ****
A+ T+ P- −0.020 0.041 −0.473 ≥0.25 ns
A+ T+ P+ −0.004 0.040 −0.104 ≥0.25 ns
Canberra_norm
(Intercept) −0.176 0.021 −8.202 <0.001 ****
A- T- P+ 0.070 0.030 2.300 0.022 *
A+ T- P- 0.040 0.031 1.307 0.192 ns
A+ T- P+ 0.083 0.031 2.726 0.007 **
A- T+ P- 0.179 0.031 5.696 <0.001 ****
A- T+ P+ 0.054 0.031 1.742 0.082 ns
A+ T+ P- 0.155 0.033 4.698 <0.001 ****
A+ T+ P+ 0.096 0.032 3.044 0.002 **
ANOVA of GLM - Beta
ANOVA(GLM(Displacement ~ Treatment), type = 2); All treatments
term statistic df p.value sig
Bray-Curtis_norm
Treatment 70.305 7.000 <0.001 ****
Canberra_norm
Treatment 46.734 7.000 <0.001 ****
Pairwise Tukey's HSD, p.adj: Dunnett - Beta
Tukey(Displacement ~ Treatment); All treatments
term .y. group1 group2 estimate std.error statistic adj.p.value Variable
Bray-Curtis_norm
Treatment Beta.Score A- T- P+ A- T- P- 0.142 0.038 3.716 0.005 Treatment
Treatment Beta.Score A+ T- P- A- T- P- 0.016 0.039 0.412 ≥0.25 Treatment
Treatment Beta.Score A+ T- P+ A- T- P- 0.002 0.039 0.040 ≥0.25 Treatment
Treatment Beta.Score A- T+ P- A- T- P- 0.081 0.040 2.048 ≥0.25 Treatment
Treatment Beta.Score A- T+ P+ A- T- P- −0.163 0.039 −4.216 <0.001 Treatment
Treatment Beta.Score A+ T+ P- A- T- P- −0.020 0.041 −0.473 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A- T- P- −0.004 0.040 −0.104 ≥0.25 Treatment
Treatment Beta.Score A+ T- P- A- T- P+ −0.126 0.039 −3.232 0.027 Treatment
Treatment Beta.Score A+ T- P+ A- T- P+ −0.141 0.039 −3.652 0.006 Treatment
Treatment Beta.Score A- T+ P- A- T- P+ −0.061 0.040 −1.543 ≥0.25 Treatment
Treatment Beta.Score A- T+ P+ A- T- P+ −0.306 0.039 −7.884 <0.001 Treatment
Treatment Beta.Score A+ T+ P- A- T- P+ −0.162 0.041 −3.900 0.002 Treatment
Treatment Beta.Score A+ T+ P+ A- T- P+ −0.146 0.040 −3.682 0.006 Treatment
Treatment Beta.Score A+ T- P+ A+ T- P- −0.015 0.039 −0.370 ≥0.25 Treatment
Treatment Beta.Score A- T+ P- A+ T- P- 0.065 0.040 1.612 ≥0.25 Treatment
Treatment Beta.Score A- T+ P+ A+ T- P- −0.180 0.040 −4.543 <0.001 Treatment
Treatment Beta.Score A+ T+ P- A+ T- P- −0.036 0.042 −0.846 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A+ T- P- −0.020 0.040 −0.499 ≥0.25 Treatment
Treatment Beta.Score A- T+ P- A+ T- P+ 0.080 0.040 1.997 ≥0.25 Treatment
Treatment Beta.Score A- T+ P+ A+ T- P+ −0.165 0.039 −4.230 <0.001 Treatment
Treatment Beta.Score A+ T+ P- A+ T- P+ −0.021 0.042 −0.507 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A+ T- P+ −0.006 0.040 −0.142 ≥0.25 Treatment
Treatment Beta.Score A- T+ P+ A- T+ P- −0.245 0.040 −6.102 <0.001 Treatment
Treatment Beta.Score A+ T+ P- A- T+ P- −0.101 0.043 −2.358 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A- T+ P- −0.085 0.041 −2.077 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P- A- T+ P+ 0.144 0.042 3.429 0.014 Treatment
Treatment Beta.Score A+ T+ P+ A- T+ P+ 0.159 0.040 3.961 0.002 Treatment
Treatment Beta.Score A+ T+ P+ A+ T+ P- 0.015 0.043 0.361 ≥0.25 Treatment
Canberra_norm
Treatment Beta.Score A- T- P+ A- T- P- 0.070 0.030 2.300 ≥0.25 Treatment
Treatment Beta.Score A+ T- P- A- T- P- 0.040 0.031 1.307 ≥0.25 Treatment
Treatment Beta.Score A+ T- P+ A- T- P- 0.083 0.031 2.726 0.115 Treatment
Treatment Beta.Score A- T+ P- A- T- P- 0.179 0.031 5.696 <0.001 Treatment
Treatment Beta.Score A- T+ P+ A- T- P- 0.054 0.031 1.742 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P- A- T- P- 0.155 0.033 4.698 <0.001 Treatment
Treatment Beta.Score A+ T+ P+ A- T- P- 0.096 0.032 3.044 0.048 Treatment
Treatment Beta.Score A+ T- P- A- T- P+ −0.029 0.031 −0.948 ≥0.25 Treatment
Treatment Beta.Score A+ T- P+ A- T- P+ 0.013 0.031 0.441 ≥0.25 Treatment
Treatment Beta.Score A- T+ P- A- T- P+ 0.109 0.031 3.473 0.012 Treatment
Treatment Beta.Score A- T+ P+ A- T- P+ −0.016 0.031 −0.528 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P- A- T- P+ 0.085 0.033 2.576 0.164 Treatment
Treatment Beta.Score A+ T+ P+ A- T- P+ 0.026 0.032 0.830 ≥0.25 Treatment
Treatment Beta.Score A+ T- P+ A+ T- P- 0.043 0.031 1.374 ≥0.25 Treatment
Treatment Beta.Score A- T+ P- A+ T- P- 0.138 0.032 4.327 <0.001 Treatment
Treatment Beta.Score A- T+ P+ A+ T- P- 0.013 0.031 0.418 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P- A+ T- P- 0.114 0.033 3.410 0.015 Treatment
Treatment Beta.Score A+ T+ P+ A+ T- P- 0.055 0.032 1.729 ≥0.25 Treatment
Treatment Beta.Score A- T+ P- A+ T- P+ 0.096 0.032 3.027 0.050 Treatment
Treatment Beta.Score A- T+ P+ A+ T- P+ −0.030 0.031 −0.960 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P- A+ T- P+ 0.071 0.033 2.156 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A+ T- P+ 0.013 0.032 0.400 ≥0.25 Treatment
Treatment Beta.Score A- T+ P+ A- T+ P- −0.125 0.032 −3.943 0.002 Treatment
Treatment Beta.Score A+ T+ P- A- T+ P- −0.024 0.034 −0.718 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A- T+ P- −0.083 0.033 −2.549 0.175 Treatment
Treatment Beta.Score A+ T+ P- A- T+ P+ 0.101 0.033 3.036 0.049 Treatment
Treatment Beta.Score A+ T+ P+ A- T+ P+ 0.042 0.032 1.329 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A+ T+ P- −0.059 0.034 −1.725 ≥0.25 Treatment

60 DPE

Displacement
Plot

Code
# Set the subsection for this analysis
tmp.resSubSection <- "All_60DPE"

# GLM Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]] <- 
  test.beta.all_metrics$displacement %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models_beta(formula_str = "Displacement ~ Treatment", 
                   beta_metric_col = "Beta.Metric", 
                   beta_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["ANOVA"]] <-
  run_glm_anova_beta(beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]])

# Tukey Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["Tukey"]] <-
  test.beta.all_metrics$displacement %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm_beta(., "Displacement", "Beta.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
Tables
GLM Results - 60 DPE
glm(Displacement ~ Treatment); All treatments at 60 DPE
term estimate std.error statistic p.value p.adj.sig
Bray-Curtis_norm
(Intercept) −0.060 0.041 −1.474 0.142 ns
A- T- P+ 0.284 0.057 4.945 <0.001 ****
A+ T- P- 0.007 0.061 0.115 ≥0.25 ns
A+ T- P+ −0.021 0.058 −0.360 ≥0.25 ns
A- T+ P- 0.043 0.063 0.675 ≥0.25 ns
A- T+ P+ −0.206 0.060 −3.458 <0.001 ***
A+ T+ P- −0.119 0.075 −1.584 0.115 ns
A+ T+ P+ 0.063 0.064 0.993 ≥0.25 ns
Canberra_norm
(Intercept) −0.124 0.034 −3.647 <0.001 ***
A- T- P+ 0.134 0.048 2.796 0.006 **
A+ T- P- −0.089 0.051 −1.750 0.081 ns
A+ T- P+ 0.037 0.049 0.765 ≥0.25 ns
A- T+ P- 0.099 0.053 1.871 0.063 ns
A- T+ P+ 0.065 0.050 1.306 0.193 ns
A+ T+ P- 0.123 0.063 1.962 0.051 ns
A+ T+ P+ 0.071 0.053 1.329 0.185 ns
ANOVA of GLM - 60 DPE
ANOVA(GLM(Displacement ~ Treatment), type = 2); All treatments at 60 DPE
term statistic df p.value sig
Bray-Curtis_norm
Treatment 77.159 7.000 <0.001 ****
Canberra_norm
Treatment 26.083 7.000 <0.001 ***
Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); All treatments at 60 DPE
term .y. group1 group2 estimate std.error statistic adj.p.value Variable
Bray-Curtis_norm
Treatment Beta.Score A- T- P+ A- T- P- 0.284 0.057 4.945 <0.001 Treatment
Treatment Beta.Score A+ T- P- A- T- P- 0.007 0.061 0.115 ≥0.25 Treatment
Treatment Beta.Score A+ T- P+ A- T- P- −0.021 0.058 −0.360 ≥0.25 Treatment
Treatment Beta.Score A- T+ P- A- T- P- 0.043 0.063 0.675 ≥0.25 Treatment
Treatment Beta.Score A- T+ P+ A- T- P- −0.206 0.060 −3.458 0.013 Treatment
Treatment Beta.Score A+ T+ P- A- T- P- −0.119 0.075 −1.584 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A- T- P- 0.063 0.064 0.993 ≥0.25 Treatment
Treatment Beta.Score A+ T- P- A- T- P+ −0.277 0.060 −4.595 <0.001 Treatment
Treatment Beta.Score A+ T- P+ A- T- P+ −0.305 0.058 −5.272 <0.001 Treatment
Treatment Beta.Score A- T+ P- A- T- P+ −0.241 0.063 −3.845 0.003 Treatment
Treatment Beta.Score A- T+ P+ A- T- P+ −0.490 0.059 −8.275 <0.001 Treatment
Treatment Beta.Score A+ T+ P- A- T- P+ −0.403 0.075 −5.372 <0.001 Treatment
Treatment Beta.Score A+ T+ P+ A- T- P+ −0.220 0.063 −3.473 0.012 Treatment
Treatment Beta.Score A+ T- P+ A+ T- P- −0.028 0.061 −0.458 ≥0.25 Treatment
Treatment Beta.Score A- T+ P- A+ T- P- 0.036 0.066 0.542 ≥0.25 Treatment
Treatment Beta.Score A- T+ P+ A+ T- P- −0.213 0.062 −3.418 0.014 Treatment
Treatment Beta.Score A+ T+ P- A+ T- P- −0.126 0.078 −1.629 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A+ T- P- 0.056 0.066 0.849 ≥0.25 Treatment
Treatment Beta.Score A- T+ P- A+ T- P+ 0.064 0.063 1.001 ≥0.25 Treatment
Treatment Beta.Score A- T+ P+ A+ T- P+ −0.185 0.060 −3.086 0.042 Treatment
Treatment Beta.Score A+ T+ P- A+ T- P+ −0.098 0.076 −1.300 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A+ T- P+ 0.084 0.064 1.313 ≥0.25 Treatment
Treatment Beta.Score A- T+ P+ A- T+ P- −0.249 0.065 −3.840 0.003 Treatment
Treatment Beta.Score A+ T+ P- A- T+ P- −0.162 0.079 −2.037 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A- T+ P- 0.021 0.069 0.303 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P- A- T+ P+ 0.087 0.077 1.130 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A- T+ P+ 0.269 0.065 4.116 <0.001 Treatment
Treatment Beta.Score A+ T+ P+ A+ T+ P- 0.183 0.080 2.282 ≥0.25 Treatment
Canberra_norm
Treatment Beta.Score A- T- P+ A- T- P- 0.134 0.048 2.796 0.094 Treatment
Treatment Beta.Score A+ T- P- A- T- P- −0.089 0.051 −1.750 ≥0.25 Treatment
Treatment Beta.Score A+ T- P+ A- T- P- 0.037 0.049 0.765 ≥0.25 Treatment
Treatment Beta.Score A- T+ P- A- T- P- 0.099 0.053 1.871 ≥0.25 Treatment
Treatment Beta.Score A- T+ P+ A- T- P- 0.065 0.050 1.306 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P- A- T- P- 0.123 0.063 1.962 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A- T- P- 0.071 0.053 1.329 ≥0.25 Treatment
Treatment Beta.Score A+ T- P- A- T- P+ −0.223 0.050 −4.425 <0.001 Treatment
Treatment Beta.Score A+ T- P+ A- T- P+ −0.097 0.048 −2.006 ≥0.25 Treatment
Treatment Beta.Score A- T+ P- A- T- P+ −0.035 0.052 −0.676 ≥0.25 Treatment
Treatment Beta.Score A- T+ P+ A- T- P+ −0.069 0.049 −1.397 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P- A- T- P+ −0.011 0.063 −0.168 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A- T- P+ −0.063 0.053 −1.192 ≥0.25 Treatment
Treatment Beta.Score A+ T- P+ A+ T- P- 0.126 0.051 2.468 0.206 Treatment
Treatment Beta.Score A- T+ P- A+ T- P- 0.187 0.055 3.412 0.014 Treatment
Treatment Beta.Score A- T+ P+ A+ T- P- 0.153 0.052 2.951 0.062 Treatment
Treatment Beta.Score A+ T+ P- A+ T- P- 0.212 0.065 3.275 0.023 Treatment
Treatment Beta.Score A+ T+ P+ A+ T- P- 0.159 0.055 2.875 0.076 Treatment
Treatment Beta.Score A- T+ P- A+ T- P+ 0.061 0.053 1.159 ≥0.25 Treatment
Treatment Beta.Score A- T+ P+ A+ T- P+ 0.028 0.050 0.555 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P- A+ T- P+ 0.086 0.063 1.366 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A+ T- P+ 0.034 0.054 0.628 ≥0.25 Treatment
Treatment Beta.Score A- T+ P+ A- T+ P- −0.034 0.054 −0.623 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P- A- T+ P- 0.025 0.066 0.375 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A- T+ P- −0.028 0.057 −0.485 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P- A- T+ P+ 0.058 0.064 0.913 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A- T+ P+ 0.006 0.055 0.107 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A+ T+ P- −0.053 0.067 −0.788 ≥0.25 Treatment
Significant Results:Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); All treatments at 60 DPE
term .y. group1 group2 estimate std.error statistic adj.p.value Variable
Bray-Curtis_norm
Treatment Beta.Score A- T- P+ A- T- P- 0.284 0.057 4.945 <0.001 Treatment
Treatment Beta.Score A- T+ P+ A- T- P- −0.206 0.060 −3.458 0.013 Treatment
Treatment Beta.Score A+ T- P- A- T- P+ −0.277 0.060 −4.595 <0.001 Treatment
Treatment Beta.Score A+ T- P+ A- T- P+ −0.305 0.058 −5.272 <0.001 Treatment
Treatment Beta.Score A- T+ P- A- T- P+ −0.241 0.063 −3.845 0.003 Treatment
Treatment Beta.Score A- T+ P+ A- T- P+ −0.490 0.059 −8.275 <0.001 Treatment
Treatment Beta.Score A+ T+ P- A- T- P+ −0.403 0.075 −5.372 <0.001 Treatment
Treatment Beta.Score A+ T+ P+ A- T- P+ −0.220 0.063 −3.473 0.012 Treatment
Treatment Beta.Score A- T+ P+ A+ T- P- −0.213 0.062 −3.418 0.014 Treatment
Treatment Beta.Score A- T+ P+ A+ T- P+ −0.185 0.060 −3.086 0.042 Treatment
Treatment Beta.Score A- T+ P+ A- T+ P- −0.249 0.065 −3.840 0.003 Treatment
Treatment Beta.Score A+ T+ P+ A- T+ P+ 0.269 0.065 4.116 <0.001 Treatment
Canberra_norm
Treatment Beta.Score A+ T- P- A- T- P+ −0.223 0.050 −4.425 <0.001 Treatment
Treatment Beta.Score A- T+ P- A+ T- P- 0.187 0.055 3.412 0.014 Treatment
Treatment Beta.Score A+ T+ P- A+ T- P- 0.212 0.065 3.275 0.023 Treatment
PCoA
Plot

Code
# Set the subsection for this analysis
tmp.resSubSection <- "All_60DPE"



#### Dispersion --------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]] <-
  run_BetaDispersion(dist.matrix = beta.dist.mat[[tmp.resSubSection]], 
                     var = c("treatment_group"))

##### ANOVA ------------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA"]] <-
  get_HoD_anova(betaDisper = beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]],
                var = c("treatment_group"))

### Table

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA.Table"]] <-
  beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA"]] %>%
  
    # Create the table
    dplyr::group_by(Beta.Metric) %>%
    set_GT(var = "p.value", group.by = "Beta.Metric")  %>%
  
    # Title/caption
    gt::tab_header(
      title = "ANOVA: Homogeneity of Dispersion",
      subtitle = "ANOVA(Beta Disperson ~ Treatment); All treatments at 60 DPE"
    ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = c("p.value", "sig"),
      rows = p.value < 0.05
    )
  )

##### Tukey ------------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey"]] <-
  get_HoD_tukey(betaDisper = beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]],
                var = c("treatment_group")) 

### Table

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey.Table"]] <-
  beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey"]] %>%
    
    dplyr::mutate(
        group1 = case_when(
        group1 == "Control" ~ "A- T- P-",
        group1 == "Parasite" ~ "A- T- P+",
        group1 == "Antibiotics" ~ "A+ T- P-",
        group1 == "Antibiotics_Parasite" ~ "A+ T- P+",
        group1 == "Temperature" ~ "A- T+ P-",
        group1 == "Temperature_Parasite" ~ "A- T+ P+",
        group1 == "Antibiotics_Temperature" ~ "A+ T+ P-",
        group1 == "Antibiotics_Temperature_Parasite" ~ "A+ T+ P+",
    ),
        group2 = case_when(
        group2 == "Control" ~ "A- T- P-",
        group2 == "Parasite" ~ "A- T- P+",
        group2 == "Antibiotics" ~ "A+ T- P-",
        group2 == "Antibiotics_Parasite" ~ "A+ T- P+",
        group2 == "Temperature" ~ "A- T+ P-",
        group2 == "Temperature_Parasite" ~ "A- T+ P+",
        group2 == "Antibiotics_Temperature" ~ "A+ T+ P-",
        group2 == "Antibiotics_Temperature_Parasite" ~ "A+ T+ P+",
    )) %>%
  
  # Create the table
  dplyr::group_by(Beta.Metric) %>%
  set_GT(var = "adj.p.value", group.by = "Beta.Metric")  %>%
  
  # Title/caption
  gt::tab_header(
    title = "Tukey: Homogeneity of Dispersion",
    subtitle = "Tukey(Beta Disperson ~ Treatment); All treatments at 60 DPE"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = c("adj.p.value", "sig"),
      rows = adj.p.value < 0.05
    )
  )
Tables
ANOVA: Homogeneity of Dispersion
ANOVA(Beta Disperson ~ Treatment); All treatments at 60 DPE
term df sumsq meansq statistic p.value sig
bray
treatment_group 7.000 0.645 0.092 3.523 0.001 **
Residual 228.000 5.965 0.026 NA NA NA
canberra
treatment_group 7.000 0.132 0.019 5.021 <0.001 ****
Residual 228.000 0.857 0.004 NA NA NA
Tukey: Homogeneity of Dispersion
Tukey(Beta Disperson ~ Treatment); All treatments at 60 DPE
.y. term group1 group2 estimate conf.low conf.high adj.p.value sig
bray
Distance treatment_group A- T- P+ A- T- P- −0.025 −0.141 0.091 ≥0.25 ns
Distance treatment_group A+ T- P- A- T- P- −0.036 −0.158 0.086 ≥0.25 ns
Distance treatment_group A+ T- P+ A- T- P- −0.065 −0.183 0.052 ≥0.25 ns
Distance treatment_group A- T+ P- A- T- P- −0.064 −0.191 0.064 ≥0.25 ns
Distance treatment_group A- T+ P+ A- T- P- −0.106 −0.227 0.014 0.126 ns
Distance treatment_group A+ T+ P- A- T- P- −0.162 −0.314 −0.010 0.028 *
Distance treatment_group A+ T+ P+ A- T- P- −0.157 −0.286 −0.028 0.006 **
Distance treatment_group A+ T- P- A- T- P+ −0.011 −0.132 0.111 ≥0.25 ns
Distance treatment_group A+ T- P+ A- T- P+ −0.040 −0.157 0.076 ≥0.25 ns
Distance treatment_group A- T+ P- A- T- P+ −0.039 −0.165 0.088 ≥0.25 ns
Distance treatment_group A- T+ P+ A- T- P+ −0.081 −0.201 0.038 ≥0.25 ns
Distance treatment_group A+ T+ P- A- T- P+ −0.137 −0.289 0.014 0.108 ns
Distance treatment_group A+ T+ P+ A- T- P+ −0.132 −0.260 −0.004 0.038 *
Distance treatment_group A+ T- P+ A+ T- P- −0.029 −0.152 0.094 ≥0.25 ns
Distance treatment_group A- T+ P- A+ T- P- −0.028 −0.160 0.105 ≥0.25 ns
Distance treatment_group A- T+ P+ A+ T- P- −0.070 −0.196 0.055 ≥0.25 ns
Distance treatment_group A+ T+ P- A+ T- P- −0.126 −0.283 0.030 0.215 ns
Distance treatment_group A+ T+ P+ A+ T- P- −0.121 −0.255 0.013 0.108 ns
Distance treatment_group A- T+ P- A+ T- P+ 0.001 −0.127 0.130 ≥0.25 ns
Distance treatment_group A- T+ P+ A+ T- P+ −0.041 −0.162 0.080 ≥0.25 ns
Distance treatment_group A+ T+ P- A+ T- P+ −0.097 −0.250 0.056 ≥0.25 ns
Distance treatment_group A+ T+ P+ A+ T- P+ −0.092 −0.222 0.038 ≥0.25 ns
Distance treatment_group A- T+ P+ A- T+ P- −0.042 −0.173 0.088 ≥0.25 ns
Distance treatment_group A+ T+ P- A- T+ P- −0.098 −0.259 0.062 ≥0.25 ns
Distance treatment_group A+ T+ P+ A- T+ P- −0.093 −0.232 0.045 ≥0.25 ns
Distance treatment_group A+ T+ P- A- T+ P+ −0.056 −0.211 0.099 ≥0.25 ns
Distance treatment_group A+ T+ P+ A- T+ P+ −0.051 −0.183 0.081 ≥0.25 ns
Distance treatment_group A+ T+ P+ A+ T+ P- 0.005 −0.157 0.166 ≥0.25 ns
canberra
Distance treatment_group A- T- P+ A- T- P- 0.008 −0.036 0.052 ≥0.25 ns
Distance treatment_group A+ T- P- A- T- P- −0.056 −0.103 −0.010 0.006 **
Distance treatment_group A+ T- P+ A- T- P- 0.007 −0.037 0.052 ≥0.25 ns
Distance treatment_group A- T+ P- A- T- P- 0.010 −0.039 0.058 ≥0.25 ns
Distance treatment_group A- T+ P+ A- T- P- −0.007 −0.052 0.039 ≥0.25 ns
Distance treatment_group A+ T+ P- A- T- P- −0.023 −0.080 0.035 ≥0.25 ns
Distance treatment_group A+ T+ P+ A- T- P- −0.044 −0.093 0.005 0.116 ns
Distance treatment_group A+ T- P- A- T- P+ −0.065 −0.111 −0.019 <0.001 ***
Distance treatment_group A+ T- P+ A- T- P+ −0.001 −0.045 0.043 ≥0.25 ns
Distance treatment_group A- T+ P- A- T- P+ 0.001 −0.047 0.049 ≥0.25 ns
Distance treatment_group A- T+ P+ A- T- P+ −0.015 −0.060 0.030 ≥0.25 ns
Distance treatment_group A+ T+ P- A- T- P+ −0.031 −0.088 0.026 ≥0.25 ns
Distance treatment_group A+ T+ P+ A- T- P+ −0.052 −0.101 −0.004 0.026 *
Distance treatment_group A+ T- P+ A+ T- P- 0.064 0.017 0.110 0.001 **
Distance treatment_group A- T+ P- A+ T- P- 0.066 0.016 0.116 0.002 **
Distance treatment_group A- T+ P+ A+ T- P- 0.050 0.002 0.097 0.034 *
Distance treatment_group A+ T+ P- A+ T- P- 0.034 −0.026 0.093 ≥0.25 ns
Distance treatment_group A+ T+ P+ A+ T- P- 0.013 −0.038 0.063 ≥0.25 ns
Distance treatment_group A- T+ P- A+ T- P+ 0.002 −0.046 0.051 ≥0.25 ns
Distance treatment_group A- T+ P+ A+ T- P+ −0.014 −0.060 0.032 ≥0.25 ns
Distance treatment_group A+ T+ P- A+ T- P+ −0.030 −0.088 0.028 ≥0.25 ns
Distance treatment_group A+ T+ P+ A+ T- P+ −0.051 −0.100 −0.002 0.035 *
Distance treatment_group A- T+ P+ A- T+ P- −0.016 −0.066 0.033 ≥0.25 ns
Distance treatment_group A+ T+ P- A- T+ P- −0.032 −0.093 0.029 ≥0.25 ns
Distance treatment_group A+ T+ P+ A- T+ P- −0.053 −0.106 −0.001 0.044 *
Distance treatment_group A+ T+ P- A- T+ P+ −0.016 −0.075 0.043 ≥0.25 ns
Distance treatment_group A+ T+ P+ A- T+ P+ −0.037 −0.087 0.013 ≥0.25 ns
Distance treatment_group A+ T+ P+ A+ T+ P- −0.021 −0.082 0.040 ≥0.25 ns
CAP
Plot

Code
# Set the subsection for this analysis
tmp.resSubSection <- "All_60DPE"


#### Capscale ----------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.mod"]] <-
  run_capscale(ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          # microViz::ps_filter(Treatment  %in% c(
          # "A- T- P-",
          # "A- T- P+",
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          # )) %>%
          microViz::ps_filter(Time == 60), 
               dist.matrix = beta.dist.mat[[tmp.resSubSection]], 
               formula_str = "dist ~ Treatment")

##### ADONIS ------------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS"]] <-
  run_cap_adonis(ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          # microViz::ps_filter(Treatment  %in% c(
          # "A- T- P-",
          # "A- T- P+",
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          # )) %>%
          microViz::ps_filter(Time == 60),
                 dist.matrix = beta.dist.mat[[tmp.resSubSection]], 
                 formula_str = "dist ~ Treatment",
                 by.method = "terms") 
## Warning in tidy.anova(.): The following column names in ANOVA output were not
## recognized or transformed: SumOfSqs, R2

## Warning in tidy.anova(.): The following column names in ANOVA output were not
## recognized or transformed: SumOfSqs, R2
### Table

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS.Table"]] <-
  beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS"]] %>%
    set_GT(var = "p.value", group.by = "Beta.Metric") %>%
  
  # Title/caption
  gt::tab_header(
    title = "ADONIS2",
    subtitle = "adonis2(Beta Distance ~ Treatment); All treatments at 60 DPE"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = c("p.value", "sig"),
      rows = p.value < 0.05
    )
  )
Tables
ADONIS2
adonis2(Beta Distance ~ Treatment); All treatments at 60 DPE
term df SumOfSqs R2 statistic p.value sig
bray
Treatment 7.000 3.508 0.105 3.813 0.001 ***
Residual 228.000 29.963 0.895 NA NA NA
Total 235.000 33.471 1.000 NA NA NA
canberra
Treatment 7.000 8.648 0.103 3.733 0.001 ***
Residual 228.000 75.464 0.897 NA NA NA
Total 235.000 84.112 1.000 NA NA NA

01) What is the host-microbiome response to parasite exposure?

Alpha

All

Displacement
Plot

Code
# Set the subsection for this analysis
tmp.resSubSection <- "ParasiteExp"

# GLM Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]] <- 
  test.alpha.all_metrics$volatility %>%
    # Filter 
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          "A- T- P+"#,
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          )) %>%
    
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models(formula_str = "Displacement ~ Treatment", 
                   alpha_metric_col = "Alpha.Metric", 
                   alpha_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["ANOVA"]] <-
  run_glm_anova(alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]])

# Tukey Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["Tukey"]] <-
  test.alpha.all_metrics$volatility %>%
    # Filter 
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          "A- T- P+"#,
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          )) %>%
    
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm(., "Displacement", "Alpha.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 
Tables
GLM Results
glm(Displacement ~ Treatment); Parasite exposure only
term estimate std.error statistic p.value p.adj.sig
Shannon
(Intercept) 0.078 0.021 3.755 <0.001 ***
A- T- P+ −0.045 0.029 −1.529 0.128 ns
Simpson
(Intercept) 0.042 0.029 1.480 0.141 ns
A- T- P+ −0.049 0.040 −1.199 0.232 ns
Richness
(Intercept) 0.221 0.025 8.745 <0.001 ****
A- T- P+ 0.009 0.036 0.251 ≥0.25 ns
ANOVA of GLM
ANOVA(GLM(Displacement ~ Treatment), type = 2); Parasite exposure only
term statistic df p.value sig
Shannon
Treatment 2.338 1.000 0.126 ns
Simpson
Treatment 1.439 1.000 0.230 ns
Richness
Treatment 0.063 1.000 ≥0.25 ns
Pairwise Tukey's HSD, p.adj: Dunnett
Tukey(Displacement ~ Treatment); Parasite exposure only
term .y. group1 group2 estimate std.error statistic adj.p.value Variable
Shannon
Treatment Alpha.Score A- T- P+ A- T- P- −0.045 0.029 −1.529 0.126 Treatment
Simpson
Treatment Alpha.Score A- T- P+ A- T- P- −0.049 0.040 −1.199 0.230 Treatment
Richness
Treatment Alpha.Score A- T- P+ A- T- P- 0.009 0.036 0.251 ≥0.25 Treatment

60 DPE

Displacement
Plot

Code
# Set the subsection for this analysis
tmp.resSubSection <- "ParasiteExp_60DPE"

# GLM Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]] <- 
  test.alpha.all_metrics$displacement %>%
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          "A- T- P+"#,
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          )) %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models(formula_str = "Displacement ~ Treatment", 
                   alpha_metric_col = "Alpha.Metric", 
                   alpha_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["ANOVA"]] <-
  run_glm_anova(alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]])

# Tukey Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["Tukey"]] <-
  test.alpha.all_metrics$displacement %>%
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          "A- T- P+"#,
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          )) %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm(., "Displacement", "Alpha.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 
Tables
GLM Results - 60 DPE
glm(Displacement ~ Treatment); All treatments at 60 DPE
term estimate std.error statistic p.value p.adj.sig
Shannon
(Intercept) 0.091 0.041 2.207 0.031 *
A- T- P+ −0.082 0.058 −1.428 0.158 ns
Simpson
(Intercept) 0.008 0.053 0.150 ≥0.25 ns
A- T- P+ −0.070 0.075 −0.933 ≥0.25 ns
Richness
(Intercept) 0.338 0.034 9.916 <0.001 ****
A- T- P+ 0.025 0.048 0.518 ≥0.25 ns
ANOVA of GLM - 60 DPE
ANOVA(GLM(Displacement ~ Treatment), type = 2); Parasite exposure only at 60 DPE
term statistic df p.value sig
Shannon
Treatment 2.039 1.000 0.153 ns
Simpson
Treatment 0.871 1.000 ≥0.25 ns
Richness
Treatment 0.269 1.000 ≥0.25 ns
Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); Parasite exposure only at 60 DPE
term .y. group1 group2 estimate std.error statistic adj.p.value Variable
Shannon
Treatment Alpha.Score A- T- P+ A- T- P- −0.082 0.058 −1.428 0.153 Treatment
Simpson
Treatment Alpha.Score A- T- P+ A- T- P- −0.070 0.075 −0.933 ≥0.25 Treatment
Richness
Treatment Alpha.Score A- T- P+ A- T- P- 0.025 0.048 0.518 ≥0.25 Treatment
Significant Results:Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); All treatments at 60 DPE
term .y. group1 group2 estimate std.error statistic adj.p.value Variable

Beta

All

Displacement
Plot

Code
# Set the subsection for this analysis
tmp.resSubSection <- "ParasiteExp"

# GLM Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]] <- 
  test.beta.all_metrics$volatility %>%
        # Filter 
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          "A- T- P+"#,
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          )) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models_beta(formula_str = "Displacement ~ Treatment", 
                   beta_metric_col = "Beta.Metric", 
                   beta_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["ANOVA"]] <-
  run_glm_anova_beta(beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]])

# Tukey Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["Tukey"]] <-
  test.beta.all_metrics$volatility %>%
        # Filter 
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          "A- T- P+"#,
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          )) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm_beta(., "Displacement", "Beta.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 
Tables
GLM Results - Beta
glm(Displacement ~ Treatment); Parasite exposure only
term estimate std.error statistic p.value p.adj.sig
Bray-Curtis_norm
(Intercept) −0.061 0.026 −2.369 0.019 *
A- T- P+ 0.142 0.037 3.884 <0.001 ***
Canberra_norm
(Intercept) −0.176 0.023 −7.710 <0.001 ****
A- T- P+ 0.070 0.032 2.162 0.032 *
ANOVA of GLM - Beta
ANOVA(GLM(Displacement ~ Treatment), type = 2); Parasite exposure only
term statistic df p.value sig
Bray-Curtis_norm
Treatment 15.083 1.000 <0.001 ***
Canberra_norm
Treatment 4.673 1.000 0.031 *
Pairwise Tukey's HSD, p.adj: Dunnett - Beta
Tukey(Displacement ~ Treatment); Parasite exposure only
term .y. group1 group2 estimate std.error statistic adj.p.value Variable
Bray-Curtis_norm
Treatment Beta.Score A- T- P+ A- T- P- 0.142 0.037 3.884 <0.001 Treatment
Canberra_norm
Treatment Beta.Score A- T- P+ A- T- P- 0.070 0.032 2.162 0.031 Treatment

60 DPE

Displacement
Plot

Code
# Set the subsection for this analysis
tmp.resSubSection <- "ParasiteExp_60DPE"

# GLM Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]] <- 
  test.beta.all_metrics$displacement %>%
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          "A- T- P+"#,
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          )) %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models_beta(formula_str = "Displacement ~ Treatment", 
                   beta_metric_col = "Beta.Metric", 
                   beta_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["ANOVA"]] <-
  run_glm_anova_beta(beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]])

# Tukey Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["Tukey"]] <-
  test.beta.all_metrics$displacement %>%
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          "A- T- P+"#,
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          )) %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm_beta(., "Displacement", "Beta.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 
Tables
GLM Results - 60 DPE
glm(Displacement ~ Treatment); Parasite exposure only at 60 DPE
term estimate std.error statistic p.value p.adj.sig
Bray-Curtis_norm
(Intercept) −0.060 0.036 −1.654 0.103 ns
A- T- P+ 0.284 0.051 5.548 <0.001 ****
Canberra_norm
(Intercept) −0.124 0.028 −4.457 <0.001 ****
A- T- P+ 0.134 0.039 3.417 0.001 **
ANOVA of GLM - 60 DPE
ANOVA(GLM(Displacement ~ Treatment), type = 2); Parasite exposure only at 60 DPE
term statistic df p.value sig
Bray-Curtis_norm
Treatment 30.779 1.000 <0.001 ****
Canberra_norm
Treatment 11.677 1.000 <0.001 ***
Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); Parasite exposure only at 60 DPE
term .y. group1 group2 estimate std.error statistic adj.p.value Variable
Bray-Curtis_norm
Treatment Beta.Score A- T- P+ A- T- P- 0.284 0.051 5.548 <0.001 Treatment
Canberra_norm
Treatment Beta.Score A- T- P+ A- T- P- 0.134 0.039 3.417 <0.001 Treatment
Significant Results:Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); All treatments at 60 DPE
term .y. group1 group2 estimate std.error statistic adj.p.value Variable
Bray-Curtis_norm
Treatment Beta.Score A- T- P+ A- T- P- 0.284 0.051 5.548 <0.001 Treatment
Canberra_norm
Treatment Beta.Score A- T- P+ A- T- P- 0.134 0.039 3.417 <0.001 Treatment
PCoA
Plot

Code
# Set the subsection for this analysis
tmp.resSubSection <- "ParasiteExp_60DPE"



#### Dispersion --------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]] <-
  run_BetaDispersion(dist.matrix = beta.dist.mat[[tmp.resSubSection]], 
                     var = c("treatment_group"))

##### ANOVA ------------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA"]] <-
  get_HoD_anova(betaDisper = beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]],
                var = c("treatment_group"))

### Table

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA.Table"]] <-
  beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA"]] %>%
  
    # Create the table
    dplyr::group_by(Beta.Metric) %>%
    set_GT(var = "p.value", group.by = "Beta.Metric")  %>%
  
    # Title/caption
    gt::tab_header(
      title = "ANOVA: Homogeneity of Dispersion",
      subtitle = "ANOVA(Beta Disperson ~ Treatment); Parasite exposure only at 60 DPE"
    ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = c("p.value", "sig"),
      rows = p.value < 0.05
    )
  )

##### Tukey ------------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey"]] <-
  get_HoD_tukey(betaDisper = beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]],
                var = c("treatment_group")) 

### Table

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey.Table"]] <-
  beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey"]] %>%
    
    dplyr::mutate(
        group1 = case_when(
        group1 == "Control" ~ "A- T- P-",
        group1 == "Parasite" ~ "A- T- P+",
        group1 == "Antibiotics" ~ "A+ T- P-",
        group1 == "Antibiotics_Parasite" ~ "A+ T- P+",
        group1 == "Temperature" ~ "A- T+ P-",
        group1 == "Temperature_Parasite" ~ "A- T+ P+",
        group1 == "Antibiotics_Temperature" ~ "A+ T+ P-",
        group1 == "Antibiotics_Temperature_Parasite" ~ "A+ T+ P+",
    ),
        group2 = case_when(
        group2 == "Control" ~ "A- T- P-",
        group2 == "Parasite" ~ "A- T- P+",
        group2 == "Antibiotics" ~ "A+ T- P-",
        group2 == "Antibiotics_Parasite" ~ "A+ T- P+",
        group2 == "Temperature" ~ "A- T+ P-",
        group2 == "Temperature_Parasite" ~ "A- T+ P+",
        group2 == "Antibiotics_Temperature" ~ "A+ T+ P-",
        group2 == "Antibiotics_Temperature_Parasite" ~ "A+ T+ P+",
    )) %>%
  
  # Create the table
  dplyr::group_by(Beta.Metric) %>%
  set_GT(var = "adj.p.value", group.by = "Beta.Metric")  %>%
  
  # Title/caption
  gt::tab_header(
    title = "Tukey: Homogeneity of Dispersion",
    subtitle = "Tukey(Beta Disperson ~ Treatment); Parasite exposure only 60 DPE"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = c("adj.p.value", "sig"),
      rows = adj.p.value < 0.05
    )
  )
Tables
ANOVA: Homogeneity of Dispersion
ANOVA(Beta Disperson ~ Treatment); Parasite exposure only at 60 DPE
term df sumsq meansq statistic p.value sig
bray
treatment_group 1.000 0.011 0.011 0.388 ≥0.25 ns
Residual 71.000 2.088 0.029 NA NA NA
canberra
treatment_group 1.000 0.001 0.001 0.608 ≥0.25 ns
Residual 71.000 0.150 0.002 NA NA NA
Tukey: Homogeneity of Dispersion
Tukey(Beta Disperson ~ Treatment); Parasite exposure only 60 DPE
.y. term group1 group2 estimate conf.low conf.high adj.p.value sig
bray
Distance treatment_group A- T- P+ A- T- P- −0.025 −0.105 0.055 ≥0.25 ns
canberra
Distance treatment_group A- T- P+ A- T- P- 0.008 −0.013 0.030 ≥0.25 ns
CAP
Plot

Code
# Set the subsection for this analysis
tmp.resSubSection <- "ParasiteExp_60DPE"


#### Capscale ----------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.mod"]] <-
  run_capscale(ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Treatment  %in% c(
          "A- T- P-",
          "A- T- P+"#,
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          )) %>%
          microViz::ps_filter(Time == 60), 
               dist.matrix = beta.dist.mat[[tmp.resSubSection]], 
               formula_str = "dist ~ Treatment")

##### ADONIS ------------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS"]] <-
  run_cap_adonis(ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Treatment  %in% c(
          "A- T- P-",
          "A- T- P+"#,
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          )) %>%
          microViz::ps_filter(Time == 60),
                 dist.matrix = beta.dist.mat[[tmp.resSubSection]], 
                 formula_str = "dist ~ Treatment",
                 by.method = "terms") 
## Warning in tidy.anova(.): The following column names in ANOVA output were not
## recognized or transformed: SumOfSqs, R2

## Warning in tidy.anova(.): The following column names in ANOVA output were not
## recognized or transformed: SumOfSqs, R2
### Table

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS.Table"]] <-
  beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS"]] %>%
    set_GT(var = "p.value", group.by = "Beta.Metric") %>%
  
  # Title/caption
  gt::tab_header(
    title = "ADONIS2",
    subtitle = "adonis2(Beta Distance ~ Treatment); Parasite exposure only at 60 DPE"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = c("p.value", "sig"),
      rows = p.value < 0.05
    )
  )
Tables
ADONIS2
adonis2(Beta Distance ~ Treatment); Parasite exposure only at 60 DPE
term df SumOfSqs R2 statistic p.value sig
bray
Treatment 1.000 0.587 0.047 3.468 0.007 **
Residual 71.000 12.011 0.953 NA NA NA
Total 72.000 12.597 1.000 NA NA NA
canberra
Treatment 1.000 1.287 0.050 3.741 0.001 ***
Residual 71.000 24.425 0.950 NA NA NA
Total 72.000 25.712 1.000 NA NA NA

02) Is the host-microbiome response to parasite exposure historically contingent?

Alpha

All

Displacement
Plot

Code
# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressParaExp"

# GLM Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]] <- 
  test.alpha.all_metrics$volatility %>%
    # Filter 
    dplyr::filter(Treatment  %in% c(
          # "A- T- P-",
          "A- T- P+",
          # "A+ T- P-",
          "A+ T- P+",
          # "A- T+ P-",
          "A- T+ P+",
          # "A+ T+ P-",
          "A+ T+ P+"
          )) %>%
    
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models(formula_str = "Displacement ~ Treatment", 
                   alpha_metric_col = "Alpha.Metric", 
                   alpha_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["ANOVA"]] <-
  run_glm_anova(alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]])

# Tukey Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["Tukey"]] <-
  test.alpha.all_metrics$volatility %>%
    # Filter 
    dplyr::filter(Treatment  %in% c(
          # "A- T- P-",
          "A- T- P+",
          # "A+ T- P-",
          "A+ T- P+",
          # "A- T+ P-",
          "A- T+ P+",
          # "A+ T+ P-",
          "A+ T+ P+"
          )) %>%
    
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm(., "Displacement", "Alpha.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 
Tables
GLM Results
glm(Displacement ~ Treatment); Prior stressors and parasite exposure
term estimate std.error statistic p.value p.adj.sig
Shannon
(Intercept) 0.033 0.020 1.674 0.095 ns
A+ T- P+ −0.027 0.028 −0.956 ≥0.25 ns
A- T+ P+ −0.047 0.028 −1.675 0.095 ns
A+ T+ P+ −0.017 0.029 −0.575 ≥0.25 ns
Simpson
(Intercept) −0.006 0.028 −0.219 ≥0.25 ns
A+ T- P+ −0.037 0.040 −0.916 ≥0.25 ns
A- T+ P+ −0.105 0.041 −2.592 0.010 *
A+ T+ P+ −0.042 0.042 −1.020 ≥0.25 ns
Richness
(Intercept) 0.230 0.025 9.206 <0.001 ****
A+ T- P+ −0.054 0.036 −1.519 0.130 ns
A- T+ P+ 0.007 0.036 0.186 ≥0.25 ns
A+ T+ P+ −0.063 0.037 −1.711 0.088 ns
ANOVA of GLM
ANOVA(GLM(Displacement ~ Treatment), type = 2); Prior stressors and parasite exposure
term statistic df p.value sig
Shannon
Treatment 2.931 3.000 ≥0.25 ns
Simpson
Treatment 6.896 3.000 0.075 ns
Richness
Treatment 5.804 3.000 0.122 ns
Pairwise Tukey's HSD, p.adj: Dunnett
Tukey(Displacement ~ Treatment); Prior stressors and parasite exposure
term .y. group1 group2 estimate std.error statistic adj.p.value Variable
Shannon
Treatment Alpha.Score A+ T- P+ A- T- P+ −0.027 0.028 −0.956 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A- T- P+ −0.047 0.028 −1.675 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A- T- P+ −0.017 0.029 −0.575 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A+ T- P+ −0.020 0.028 −0.721 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A+ T- P+ 0.010 0.029 0.349 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A- T+ P+ 0.031 0.029 1.046 ≥0.25 Treatment
Simpson
Treatment Alpha.Score A+ T- P+ A- T- P+ −0.037 0.040 −0.916 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A- T- P+ −0.105 0.041 −2.592 0.047 Treatment
Treatment Alpha.Score A+ T+ P+ A- T- P+ −0.042 0.042 −1.020 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A+ T- P+ −0.068 0.041 −1.671 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A+ T- P+ −0.005 0.042 −0.131 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A- T+ P+ 0.063 0.042 1.490 ≥0.25 Treatment
Richness
Treatment Alpha.Score A+ T- P+ A- T- P+ −0.054 0.036 −1.519 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A- T- P+ 0.007 0.036 0.186 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A- T- P+ −0.063 0.037 −1.711 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A+ T- P+ 0.061 0.036 1.685 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A+ T- P+ −0.009 0.037 −0.237 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A- T+ P+ −0.070 0.037 −1.870 0.241 Treatment

60 DPE

Displacement
Plot

Code
# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressParaExp_60DPE"

# GLM Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]] <- 
  test.alpha.all_metrics$displacement %>%
    dplyr::filter(Treatment  %in% c(
          # "A- T- P-",
          "A- T- P+",
          # "A+ T- P-",
          "A+ T- P+",
          # "A- T+ P-",
          "A- T+ P+",
          # "A+ T+ P-",
          "A+ T+ P+"
          )) %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models(formula_str = "Displacement ~ Treatment", 
                   alpha_metric_col = "Alpha.Metric", 
                   alpha_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["ANOVA"]] <-
  run_glm_anova(alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]])

# Tukey Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["Tukey"]] <-
  test.alpha.all_metrics$displacement %>%
    dplyr::filter(Treatment  %in% c(
          # "A- T- P-",
          "A- T- P+",
          # "A+ T- P-",
          "A+ T- P+",
          # "A- T+ P-",
          "A- T+ P+",
          # "A+ T+ P-",
          "A+ T+ P+"
          )) %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm(., "Displacement", "Alpha.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 
Tables
GLM Results - 60 DPE
glm(Displacement ~ Treatment); All treatments at 60 DPE
term estimate std.error statistic p.value p.adj.sig
Shannon
(Intercept) 0.008 0.037 0.224 ≥0.25 ns
A+ T- P+ −0.001 0.053 −0.022 ≥0.25 ns
A- T+ P+ −0.077 0.054 −1.424 0.157 ns
A+ T+ P+ 0.004 0.058 0.065 ≥0.25 ns
Simpson
(Intercept) −0.062 0.048 −1.286 0.201 ns
A+ T- P+ 0.000 0.069 −0.005 ≥0.25 ns
A- T+ P+ −0.180 0.071 −2.547 0.012 *
A+ T+ P+ −0.037 0.076 −0.484 ≥0.25 ns
Richness
(Intercept) 0.363 0.032 11.248 <0.001 ****
A+ T- P+ −0.132 0.046 −2.863 0.005 **
A- T+ P+ −0.003 0.047 −0.069 ≥0.25 ns
A+ T+ P+ −0.053 0.051 −1.039 ≥0.25 ns
ANOVA of GLM - 60 DPE
ANOVA(GLM(Displacement ~ Treatment), type = 2); Prior stressors and parasite exposure at 60 DPE
term statistic df p.value sig
Shannon
Treatment 2.892 3.000 ≥0.25 ns
Simpson
Treatment 8.445 3.000 0.038 *
Richness
Treatment 10.375 3.000 0.016 *
Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); Prior stressors and parasite exposure at 60 DPE
term .y. group1 group2 estimate std.error statistic adj.p.value Variable
Shannon
Treatment Alpha.Score A+ T- P+ A- T- P+ −0.001 0.053 −0.022 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A- T- P+ −0.077 0.054 −1.424 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A- T- P+ 0.004 0.058 0.065 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A+ T- P+ −0.076 0.055 −1.384 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A+ T- P+ 0.005 0.059 0.084 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A- T+ P+ 0.081 0.060 1.351 ≥0.25 Treatment
Simpson
Treatment Alpha.Score A+ T- P+ A- T- P+ 0.000 0.069 −0.005 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A- T- P+ −0.180 0.071 −2.547 0.053 Treatment
Treatment Alpha.Score A+ T+ P+ A- T- P+ −0.037 0.076 −0.484 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A+ T- P+ −0.179 0.071 −2.509 0.058 Treatment
Treatment Alpha.Score A+ T+ P+ A+ T- P+ −0.036 0.077 −0.475 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A- T+ P+ 0.143 0.078 1.834 ≥0.25 Treatment
Richness
Treatment Alpha.Score A+ T- P+ A- T- P+ −0.132 0.046 −2.863 0.022 Treatment
Treatment Alpha.Score A- T+ P+ A- T- P+ −0.003 0.047 −0.069 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A- T- P+ −0.053 0.051 −1.039 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P+ A+ T- P+ 0.129 0.048 2.692 0.036 Treatment
Treatment Alpha.Score A+ T+ P+ A+ T- P+ 0.080 0.051 1.551 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P+ A- T+ P+ −0.050 0.052 −0.946 ≥0.25 Treatment
Significant Results:Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); All treatments at 60 DPE
term .y. group1 group2 estimate std.error statistic adj.p.value Variable
Richness
Treatment Alpha.Score A+ T- P+ A- T- P+ −0.132 0.046 −2.863 0.022 Treatment
Treatment Alpha.Score A- T+ P+ A+ T- P+ 0.129 0.048 2.692 0.036 Treatment

Beta

All

Displacement
Plot

Code
# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressParaExp"

# GLM Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]] <- 
  test.beta.all_metrics$volatility %>%
        # Filter 
    dplyr::filter(Treatment  %in% c(
          # "A- T- P-",
          "A- T- P+",
          # "A+ T- P-",
          "A+ T- P+",
          # "A- T+ P-",
          "A- T+ P+",
          # "A+ T+ P-",
          "A+ T+ P+"
          )) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models_beta(formula_str = "Displacement ~ Treatment", 
                   beta_metric_col = "Beta.Metric", 
                   beta_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["ANOVA"]] <-
  run_glm_anova_beta(beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]])

# Tukey Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["Tukey"]] <-
  test.beta.all_metrics$volatility %>%
        # Filter 
    dplyr::filter(Treatment  %in% c(
          # "A- T- P-",
          "A- T- P+",
          # "A+ T- P-",
          "A+ T- P+",
          # "A- T+ P-",
          "A- T+ P+",
          # "A+ T+ P-",
          "A+ T+ P+"
          )) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm_beta(., "Displacement", "Beta.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 
Tables
GLM Results - Beta
glm(Displacement ~ Treatment); Prior stressors and parasite exposure
term estimate std.error statistic p.value p.adj.sig
Bray-Curtis_norm
(Intercept) 0.081 0.027 2.973 0.003 **
A+ T- P+ −0.141 0.039 −3.633 <0.001 ***
A- T+ P+ −0.306 0.039 −7.843 <0.001 ****
A+ T+ P+ −0.146 0.040 −3.662 <0.001 ***
Canberra_norm
(Intercept) −0.106 0.020 −5.372 <0.001 ****
A+ T- P+ 0.013 0.028 0.478 ≥0.25 ns
A- T+ P+ −0.016 0.028 −0.573 ≥0.25 ns
A+ T+ P+ 0.026 0.029 0.900 ≥0.25 ns
ANOVA of GLM - Beta
ANOVA(GLM(Displacement ~ Treatment), type = 2); Prior stressors and parasite exposure
term statistic df p.value sig
Bray-Curtis_norm
Treatment 61.580 3.000 <0.001 ****
Canberra_norm
Treatment 2.326 3.000 ≥0.25 ns
Pairwise Tukey's HSD, p.adj: Dunnett - Beta
Tukey(Displacement ~ Treatment); Prior stressors and parasite exposure
term .y. group1 group2 estimate std.error statistic adj.p.value Variable
Bray-Curtis_norm
Treatment Beta.Score A+ T- P+ A- T- P+ −0.141 0.039 −3.633 0.002 Treatment
Treatment Beta.Score A- T+ P+ A- T- P+ −0.306 0.039 −7.843 <0.001 Treatment
Treatment Beta.Score A+ T+ P+ A- T- P+ −0.146 0.040 −3.662 0.002 Treatment
Treatment Beta.Score A- T+ P+ A+ T- P+ −0.165 0.039 −4.208 <0.001 Treatment
Treatment Beta.Score A+ T+ P+ A+ T- P+ −0.006 0.040 −0.141 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A- T+ P+ 0.159 0.040 3.940 <0.001 Treatment
Canberra_norm
Treatment Beta.Score A+ T- P+ A- T- P+ 0.013 0.028 0.478 ≥0.25 Treatment
Treatment Beta.Score A- T+ P+ A- T- P+ −0.016 0.028 −0.573 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A- T- P+ 0.026 0.029 0.900 ≥0.25 Treatment
Treatment Beta.Score A- T+ P+ A+ T- P+ −0.030 0.029 −1.042 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A+ T- P+ 0.013 0.029 0.434 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A- T+ P+ 0.042 0.029 1.442 ≥0.25 Treatment

60 DPE

Displacement
Plot

Code
# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressParaExp_60DPE"

# GLM Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]] <- 
  test.beta.all_metrics$displacement %>%
    dplyr::filter(Treatment  %in% c(
          # "A- T- P-",
          "A- T- P+",
          # "A+ T- P-",
          "A+ T- P+",
          # "A- T+ P-",
          "A- T+ P+",
          # "A+ T+ P-",
          "A+ T+ P+"
          )) %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models_beta(formula_str = "Displacement ~ Treatment", 
                   beta_metric_col = "Beta.Metric", 
                   beta_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["ANOVA"]] <-
  run_glm_anova_beta(beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]])

# Tukey Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["Tukey"]] <-
  test.beta.all_metrics$displacement %>%
    dplyr::filter(Treatment  %in% c(
          # "A- T- P-",
          "A- T- P+",
          # "A+ T- P-",
          "A+ T- P+",
          # "A- T+ P-",
          "A- T+ P+",
          # "A+ T+ P-",
          "A+ T+ P+"
          )) %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm_beta(., "Displacement", "Beta.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 
Tables
GLM Results - 60 DPE
glm(Displacement ~ Treatment); Prior stressors and parasite exposure only at 60 DPE
term estimate std.error statistic p.value p.adj.sig
Bray-Curtis_norm
(Intercept) 0.224 0.038 5.823 <0.001 ****
A+ T- P+ −0.305 0.055 −5.534 <0.001 ****
A- T+ P+ −0.490 0.056 −8.688 <0.001 ****
A+ T+ P+ −0.220 0.060 −3.646 <0.001 ***
Canberra_norm
(Intercept) 0.010 0.028 0.340 ≥0.25 ns
A+ T- P+ −0.097 0.040 −2.398 0.018 *
A- T+ P+ −0.069 0.041 −1.670 0.097 ns
A+ T+ P+ −0.063 0.044 −1.425 0.157 ns
ANOVA of GLM - 60 DPE
ANOVA(GLM(Displacement ~ Treatment), type = 2); Prior stressors and parasite exposure only at 60 DPE
term statistic df p.value sig
Bray-Curtis_norm
Treatment 78.442 3.000 <0.001 ****
Canberra_norm
Treatment 6.174 3.000 0.103 ns
Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); Prior stressors and parasite exposure only at 60 DPE
term .y. group1 group2 estimate std.error statistic adj.p.value Variable
Bray-Curtis_norm
Treatment Beta.Score A+ T- P+ A- T- P+ −0.305 0.055 −5.534 <0.001 Treatment
Treatment Beta.Score A- T+ P+ A- T- P+ −0.490 0.056 −8.688 0.000 Treatment
Treatment Beta.Score A+ T+ P+ A- T- P+ −0.220 0.060 −3.646 0.002 Treatment
Treatment Beta.Score A- T+ P+ A+ T- P+ −0.185 0.057 −3.240 0.007 Treatment
Treatment Beta.Score A+ T+ P+ A+ T- P+ 0.084 0.061 1.379 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A- T+ P+ 0.269 0.062 4.321 <0.001 Treatment
Canberra_norm
Treatment Beta.Score A+ T- P+ A- T- P+ −0.097 0.040 −2.398 0.077 Treatment
Treatment Beta.Score A- T+ P+ A- T- P+ −0.069 0.041 −1.670 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A- T- P+ −0.063 0.044 −1.425 ≥0.25 Treatment
Treatment Beta.Score A- T+ P+ A+ T- P+ 0.028 0.042 0.663 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A+ T- P+ 0.034 0.045 0.750 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P+ A- T+ P+ 0.006 0.046 0.128 ≥0.25 Treatment
Significant Results:Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); All treatments at 60 DPE
term .y. group1 group2 estimate std.error statistic adj.p.value Variable
Bray-Curtis_norm
Treatment Beta.Score A+ T- P+ A- T- P+ −0.305 0.055 −5.534 <0.001 Treatment
Treatment Beta.Score A- T+ P+ A- T- P+ −0.490 0.056 −8.688 0.000 Treatment
Treatment Beta.Score A+ T+ P+ A- T- P+ −0.220 0.060 −3.646 0.002 Treatment
Treatment Beta.Score A- T+ P+ A+ T- P+ −0.185 0.057 −3.240 0.007 Treatment
Treatment Beta.Score A+ T+ P+ A- T+ P+ 0.269 0.062 4.321 <0.001 Treatment
PCoA
Plot

Code
# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressParaExp_60DPE"



#### Dispersion --------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]] <-
  run_BetaDispersion(dist.matrix = beta.dist.mat[[tmp.resSubSection]], 
                     var = c("treatment_group"))

##### ANOVA ------------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA"]] <-
  get_HoD_anova(betaDisper = beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]],
                var = c("treatment_group"))

### Table

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA.Table"]] <-
  beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA"]] %>%
  
    # Create the table
    dplyr::group_by(Beta.Metric) %>%
    set_GT(var = "p.value", group.by = "Beta.Metric")  %>%
  
    # Title/caption
    gt::tab_header(
      title = "ANOVA: Homogeneity of Dispersion",
      subtitle = "ANOVA(Beta Disperson ~ Treatment); Prior stressors and parasite exposure at 60 DPE"
    ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = c("p.value", "sig"),
      rows = p.value < 0.05
    )
  )

##### Tukey ------------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey"]] <-
  get_HoD_tukey(betaDisper = beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]],
                var = c("treatment_group")) 

### Table

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey.Table"]] <-
  beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey"]] %>%
    
    dplyr::mutate(
        group1 = case_when(
        group1 == "Control" ~ "A- T- P-",
        group1 == "Parasite" ~ "A- T- P+",
        group1 == "Antibiotics" ~ "A+ T- P-",
        group1 == "Antibiotics_Parasite" ~ "A+ T- P+",
        group1 == "Temperature" ~ "A- T+ P-",
        group1 == "Temperature_Parasite" ~ "A- T+ P+",
        group1 == "Antibiotics_Temperature" ~ "A+ T+ P-",
        group1 == "Antibiotics_Temperature_Parasite" ~ "A+ T+ P+",
    ),
        group2 = case_when(
        group2 == "Control" ~ "A- T- P-",
        group2 == "Parasite" ~ "A- T- P+",
        group2 == "Antibiotics" ~ "A+ T- P-",
        group2 == "Antibiotics_Parasite" ~ "A+ T- P+",
        group2 == "Temperature" ~ "A- T+ P-",
        group2 == "Temperature_Parasite" ~ "A- T+ P+",
        group2 == "Antibiotics_Temperature" ~ "A+ T+ P-",
        group2 == "Antibiotics_Temperature_Parasite" ~ "A+ T+ P+",
    )) %>%
  
  # Create the table
  dplyr::group_by(Beta.Metric) %>%
  set_GT(var = "adj.p.value", group.by = "Beta.Metric")  %>%
  
  # Title/caption
  gt::tab_header(
    title = "Tukey: Homogeneity of Dispersion",
    subtitle = "Tukey(Beta Disperson ~ Treatment); Prior stressors and parasite exposure at 60 DPE"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = c("adj.p.value", "sig"),
      rows = adj.p.value < 0.05
    )
  )
Tables
ANOVA: Homogeneity of Dispersion
ANOVA(Beta Disperson ~ Treatment); Prior stressors and parasite exposure at 60 DPE
term df sumsq meansq statistic p.value sig
bray
treatment_group 3.000 0.291 0.097 3.785 0.012 *
Residual 125.000 3.198 0.026 NA NA NA
canberra
treatment_group 3.000 0.049 0.016 4.464 0.005 **
Residual 125.000 0.462 0.004 NA NA NA
Tukey: Homogeneity of Dispersion
Tukey(Beta Disperson ~ Treatment); Prior stressors and parasite exposure at 60 DPE
.y. term group1 group2 estimate conf.low conf.high adj.p.value sig
bray
Distance treatment_group A+ T- P+ A- T- P+ −0.040 −0.138 0.058 ≥0.25 ns
Distance treatment_group A- T+ P+ A- T- P+ −0.081 −0.182 0.019 0.157 ns
Distance treatment_group A+ T+ P+ A- T- P+ −0.132 −0.240 −0.024 0.009 **
Distance treatment_group A- T+ P+ A+ T- P+ −0.041 −0.143 0.061 ≥0.25 ns
Distance treatment_group A+ T+ P+ A+ T- P+ −0.092 −0.201 0.017 0.130 ns
Distance treatment_group A+ T+ P+ A- T+ P+ −0.051 −0.162 0.060 ≥0.25 ns
canberra
Distance treatment_group A+ T- P+ A- T- P+ −0.001 −0.038 0.036 ≥0.25 ns
Distance treatment_group A- T+ P+ A- T- P+ −0.015 −0.053 0.023 ≥0.25 ns
Distance treatment_group A+ T+ P+ A- T- P+ −0.052 −0.093 −0.011 0.007 **
Distance treatment_group A- T+ P+ A+ T- P+ −0.014 −0.053 0.025 ≥0.25 ns
Distance treatment_group A+ T+ P+ A+ T- P+ −0.051 −0.093 −0.010 0.009 **
Distance treatment_group A+ T+ P+ A- T+ P+ −0.037 −0.079 0.005 0.106 ns
CAP
Plot

Code
# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressParaExp_60DPE"


#### Capscale ----------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.mod"]] <-
  run_capscale(ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Treatment  %in% c(
          # "A- T- P-",
          "A- T- P+",
          # "A+ T- P-",
          "A+ T- P+",
          # "A- T+ P-",
          "A- T+ P+",
          # "A+ T+ P-",
          "A+ T+ P+"
          )) %>%
          microViz::ps_filter(Time == 60), 
               dist.matrix = beta.dist.mat[[tmp.resSubSection]], 
               formula_str = "dist ~ Treatment")

##### ADONIS ------------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS"]] <-
  run_cap_adonis(ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Treatment  %in% c(
          # "A- T- P-",
          "A- T- P+",
          # "A+ T- P-",
          "A+ T- P+",
          # "A- T+ P-",
          "A- T+ P+",
          # "A+ T+ P-",
          "A+ T+ P+"
          )) %>%
          microViz::ps_filter(Time == 60),
                 dist.matrix = beta.dist.mat[[tmp.resSubSection]], 
                 formula_str = "dist ~ Treatment",
                 by.method = "terms") 
## Warning in tidy.anova(.): The following column names in ANOVA output were not
## recognized or transformed: SumOfSqs, R2

## Warning in tidy.anova(.): The following column names in ANOVA output were not
## recognized or transformed: SumOfSqs, R2
### Table

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS.Table"]] <-
  beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS"]] %>%
    set_GT(var = "p.value", group.by = "Beta.Metric") %>%
  
  # Title/caption
  gt::tab_header(
    title = "ADONIS2",
    subtitle = "adonis2(Beta Distance ~ Treatment); Prior stressors and parasite exposure at 60 DPE"
  )  %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = c("p.value", "sig"),
      rows = p.value < 0.05
    )
  )
Tables
ADONIS2
adonis2(Beta Distance ~ Treatment); Prior stressors and parasite exposure at 60 DPE
term df SumOfSqs R2 statistic p.value sig
bray
Treatment 3.000 0.708 0.045 1.959 0.037 *
Residual 125.000 15.053 0.955 NA NA NA
Total 128.000 15.761 1.000 NA NA NA
canberra
Treatment 3.000 2.852 0.064 2.833 0.001 ***
Residual 125.000 41.942 0.936 NA NA NA
Total 128.000 44.794 1.000 NA NA NA

03) Is the host-microbiome recovery historically contingent?

Alpha

All

Displacement
Plot

Code
# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressNoParaExp"

# GLM Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]] <- 
  test.alpha.all_metrics$volatility %>%
    # Filter 
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          # "A- T- P+",
          "A+ T- P-",
          # "A+ T- P+",
          "A- T+ P-",
          # "A- T+ P+",
          "A+ T+ P-"#,
          # "A+ T+ P+"
          )) %>%
    
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models(formula_str = "Displacement ~ Treatment", 
                   alpha_metric_col = "Alpha.Metric", 
                   alpha_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["ANOVA"]] <-
  run_glm_anova(alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]])

# Tukey Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["Tukey"]] <-
  test.alpha.all_metrics$volatility %>%
    # Filter 
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          # "A- T- P+",
          "A+ T- P-",
          # "A+ T- P+",
          "A- T+ P-",
          # "A- T+ P+",
          "A+ T+ P-"#,
          # "A+ T+ P+"
          )) %>%
    
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm(., "Displacement", "Alpha.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 
Tables
GLM Results
glm(Displacement ~ Treatment); Prior stressors and no parasite exposure
term estimate std.error statistic p.value p.adj.sig
Shannon
(Intercept) 0.078 0.019 4.047 <0.001 ****
A+ T- P- −0.051 0.028 −1.839 0.067 ns
A- T+ P- −0.004 0.028 −0.153 ≥0.25 ns
A+ T+ P- −0.113 0.029 −3.849 <0.001 ***
Simpson
(Intercept) 0.042 0.027 1.562 0.119 ns
A+ T- P- −0.034 0.039 −0.861 ≥0.25 ns
A- T+ P- 0.035 0.040 0.894 ≥0.25 ns
A+ T+ P- −0.181 0.042 −4.369 <0.001 ****
Richness
(Intercept) 0.221 0.023 9.748 <0.001 ****
A+ T- P- −0.049 0.033 −1.508 0.133 ns
A- T+ P- −0.016 0.033 −0.493 ≥0.25 ns
A+ T+ P- 0.029 0.035 0.825 ≥0.25 ns
ANOVA of GLM
ANOVA(GLM(Displacement ~ Treatment), type = 2); Prior stressors and no parasite exposure
term statistic df p.value sig
Shannon
Treatment 18.423 3.000 <0.001 ***
Simpson
Treatment 29.201 3.000 <0.001 ****
Richness
Treatment 5.193 3.000 0.158 ns
Pairwise Tukey's HSD, p.adj: Dunnett
Tukey(Displacement ~ Treatment); Prior stressors and no parasite exposure
term .y. group1 group2 estimate std.error statistic adj.p.value Variable
Shannon
Treatment Alpha.Score A+ T- P- A- T- P- −0.051 0.028 −1.839 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A- T- P- −0.004 0.028 −0.153 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A- T- P- −0.113 0.029 −3.849 <0.001 Treatment
Treatment Alpha.Score A- T+ P- A+ T- P- 0.047 0.029 1.629 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A+ T- P- −0.062 0.030 −2.083 0.158 Treatment
Treatment Alpha.Score A+ T+ P- A- T+ P- −0.109 0.030 −3.596 0.002 Treatment
Simpson
Treatment Alpha.Score A+ T- P- A- T- P- −0.034 0.039 −0.861 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A- T- P- 0.035 0.040 0.894 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A- T- P- −0.181 0.042 −4.369 <0.001 Treatment
Treatment Alpha.Score A- T+ P- A+ T- P- 0.069 0.040 1.711 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A+ T- P- −0.148 0.042 −3.500 0.003 Treatment
Treatment Alpha.Score A+ T+ P- A- T+ P- −0.217 0.043 −5.073 <0.001 Treatment
Richness
Treatment Alpha.Score A+ T- P- A- T- P- −0.049 0.033 −1.508 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A- T- P- −0.016 0.033 −0.493 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A- T- P- 0.029 0.035 0.825 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A+ T- P- 0.033 0.034 0.975 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A+ T- P- 0.078 0.035 2.206 0.121 Treatment
Treatment Alpha.Score A+ T+ P- A- T+ P- 0.045 0.036 1.258 ≥0.25 Treatment

60 DPE

Displacement
Plot

Code
# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressNoParaExp_60DPE"

# GLM Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]] <- 
  test.alpha.all_metrics$displacement %>%
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          # "A- T- P+",
          "A+ T- P-",
          # "A+ T- P+",
          "A- T+ P-",
          # "A- T+ P+",
          "A+ T+ P-"#,
          # "A+ T+ P+"
          )) %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models(formula_str = "Displacement ~ Treatment", 
                   alpha_metric_col = "Alpha.Metric", 
                   alpha_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["ANOVA"]] <-
  run_glm_anova(alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]])

# Tukey Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["Tukey"]] <-
  test.alpha.all_metrics$displacement %>%
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          # "A- T- P+",
          "A+ T- P-",
          # "A+ T- P+",
          "A- T+ P-",
          # "A- T+ P+",
          "A+ T+ P-"#,
          # "A+ T+ P+"
          )) %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm(., "Displacement", "Alpha.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 
Tables
GLM Results - 60 DPE
glm(Displacement ~ Treatment); Prior stressors and no parasite exposure at 60 DPE
term estimate std.error statistic p.value p.adj.sig
Shannon
(Intercept) 0.091 0.036 2.496 0.014 *
A+ T- P- −0.034 0.054 −0.627 ≥0.25 ns
A- T+ P- −0.041 0.056 −0.731 ≥0.25 ns
A+ T+ P- −0.191 0.067 −2.860 0.005 **
Simpson
(Intercept) 0.008 0.050 0.161 ≥0.25 ns
A+ T- P- 0.023 0.074 0.308 ≥0.25 ns
A- T+ P- 0.030 0.077 0.395 ≥0.25 ns
A+ T+ P- −0.332 0.092 −3.627 <0.001 ***
Richness
(Intercept) 0.338 0.031 10.867 <0.001 ****
A+ T- P- −0.064 0.046 −1.388 0.168 ns
A- T+ P- −0.089 0.048 −1.856 0.066 ns
A+ T+ P- 0.067 0.057 1.174 0.243 ns
ANOVA of GLM - 60 DPE
ANOVA(GLM(Displacement ~ Treatment), type = 2); Prior stressors and no parasite exposure at 60 DPE
term statistic df p.value sig
Shannon
Treatment 8.403 3.000 0.038 *
Simpson
Treatment 17.775 3.000 <0.001 ***
Richness
Treatment 8.617 3.000 0.035 *
Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); Prior stressors and no parasite exposure at 60 DPE
term .y. group1 group2 estimate std.error statistic adj.p.value Variable
Shannon
Treatment Alpha.Score A+ T- P- A- T- P- −0.034 0.054 −0.627 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A- T- P- −0.041 0.056 −0.731 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A- T- P- −0.191 0.067 −2.860 0.022 Treatment
Treatment Alpha.Score A- T+ P- A+ T- P- −0.007 0.058 −0.124 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A+ T- P- −0.158 0.069 −2.289 0.099 Treatment
Treatment Alpha.Score A+ T+ P- A- T+ P- −0.150 0.071 −2.131 0.142 Treatment
Simpson
Treatment Alpha.Score A+ T- P- A- T- P- 0.023 0.074 0.308 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A- T- P- 0.030 0.077 0.395 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A- T- P- −0.332 0.092 −3.627 0.002 Treatment
Treatment Alpha.Score A- T+ P- A+ T- P- 0.008 0.080 0.095 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A+ T- P- −0.355 0.094 −3.765 <0.001 Treatment
Treatment Alpha.Score A+ T+ P- A- T+ P- −0.362 0.097 −3.751 0.001 Treatment
Richness
Treatment Alpha.Score A+ T- P- A- T- P- −0.064 0.046 −1.388 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A- T- P- −0.089 0.048 −1.856 0.245 Treatment
Treatment Alpha.Score A+ T+ P- A- T- P- 0.067 0.057 1.174 ≥0.25 Treatment
Treatment Alpha.Score A- T+ P- A+ T- P- −0.025 0.050 −0.502 ≥0.25 Treatment
Treatment Alpha.Score A+ T+ P- A+ T- P- 0.131 0.059 2.226 0.115 Treatment
Treatment Alpha.Score A+ T+ P- A- T+ P- 0.156 0.061 2.586 0.047 Treatment
Significant Results:Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); All treatments at 60 DPE
term .y. group1 group2 estimate std.error statistic adj.p.value Variable
Shannon
Treatment Alpha.Score A+ T+ P- A- T- P- −0.191 0.067 −2.860 0.022 Treatment
Simpson
Treatment Alpha.Score A+ T+ P- A- T- P- −0.332 0.092 −3.627 0.002 Treatment
Treatment Alpha.Score A+ T+ P- A+ T- P- −0.355 0.094 −3.765 <0.001 Treatment
Treatment Alpha.Score A+ T+ P- A- T+ P- −0.362 0.097 −3.751 0.001 Treatment
Richness
Treatment Alpha.Score A+ T+ P- A- T+ P- 0.156 0.061 2.586 0.047 Treatment

Beta

All

Displacement
Plot

Code
# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressNoParaExp"

# GLM Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]] <- 
  test.beta.all_metrics$volatility %>%
        # Filter 
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          # "A- T- P+",
          "A+ T- P-",
          # "A+ T- P+",
          "A- T+ P-",
          # "A- T+ P+",
          "A+ T+ P-"#,
          # "A+ T+ P+"
          )) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models_beta(formula_str = "Displacement ~ Treatment", 
                   beta_metric_col = "Beta.Metric", 
                   beta_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["ANOVA"]] <-
  run_glm_anova_beta(beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]])

# Tukey Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["Tukey"]] <-
  test.beta.all_metrics$volatility %>%
        # Filter 
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          # "A- T- P+",
          "A+ T- P-",
          # "A+ T- P+",
          "A- T+ P-",
          # "A- T+ P+",
          "A+ T+ P-"#,
          # "A+ T+ P+"
          )) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm_beta(., "Displacement", "Beta.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 
Tables
GLM Results - Beta
glm(Displacement ~ Treatment); Prior stressors and no parasite exposure
term estimate std.error statistic p.value p.adj.sig
Bray-Curtis_norm
(Intercept) −0.061 0.027 −2.279 0.023 *
A+ T- P- 0.016 0.039 0.414 ≥0.25 ns
A- T+ P- 0.081 0.039 2.059 0.040 *
A+ T+ P- −0.020 0.041 −0.475 ≥0.25 ns
Canberra_norm
(Intercept) −0.176 0.023 −7.610 <0.001 ****
A+ T- P- 0.040 0.033 1.213 0.226 ns
A- T+ P- 0.179 0.034 5.285 <0.001 ****
A+ T+ P- 0.155 0.035 4.359 <0.001 ****
ANOVA of GLM - Beta
ANOVA(GLM(Displacement ~ Treatment), type = 2); Prior stressors and no parasite exposure
term statistic df p.value sig
Bray-Curtis_norm
Treatment 6.694 3.000 0.082 ns
Canberra_norm
Treatment 38.039 3.000 <0.001 ****
Pairwise Tukey's HSD, p.adj: Dunnett - Beta
Tukey(Displacement ~ Treatment); Prior stressors and no parasite exposure
term .y. group1 group2 estimate std.error statistic adj.p.value Variable
Bray-Curtis_norm
Treatment Beta.Score A+ T- P- A- T- P- 0.016 0.039 0.414 ≥0.25 Treatment
Treatment Beta.Score A- T+ P- A- T- P- 0.081 0.039 2.059 0.166 Treatment
Treatment Beta.Score A+ T+ P- A- T- P- −0.020 0.041 −0.475 ≥0.25 Treatment
Treatment Beta.Score A- T+ P- A+ T- P- 0.065 0.040 1.621 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P- A+ T- P- −0.036 0.042 −0.851 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P- A- T+ P- −0.101 0.042 −2.371 0.083 Treatment
Canberra_norm
Treatment Beta.Score A+ T- P- A- T- P- 0.040 0.033 1.213 ≥0.25 Treatment
Treatment Beta.Score A- T+ P- A- T- P- 0.179 0.034 5.285 <0.001 Treatment
Treatment Beta.Score A+ T+ P- A- T- P- 0.155 0.035 4.359 <0.001 Treatment
Treatment Beta.Score A- T+ P- A+ T- P- 0.138 0.034 4.015 <0.001 Treatment
Treatment Beta.Score A+ T+ P- A+ T- P- 0.114 0.036 3.164 0.009 Treatment
Treatment Beta.Score A+ T+ P- A- T+ P- −0.024 0.037 −0.666 ≥0.25 Treatment

60 DPE

Displacement
Plot

Code
# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressNoParaExp_60DPE"

# GLM Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]] <- 
  test.beta.all_metrics$displacement %>%
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          # "A- T- P+",
          "A+ T- P-",
          # "A+ T- P+",
          "A- T+ P-",
          # "A- T+ P+",
          "A+ T+ P-"#,
          # "A+ T+ P+"
          )) %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models_beta(formula_str = "Displacement ~ Treatment", 
                   beta_metric_col = "Beta.Metric", 
                   beta_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["ANOVA"]] <-
  run_glm_anova_beta(beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]])

# Tukey Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["Tukey"]] <-
  test.beta.all_metrics$displacement %>%
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          # "A- T- P+",
          "A+ T- P-",
          # "A+ T- P+",
          "A- T+ P-",
          # "A- T+ P+",
          "A+ T+ P-"#,
          # "A+ T+ P+"
          )) %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm_beta(., "Displacement", "Beta.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 
Tables
GLM Results - 60 DPE
glm(Displacement ~ Treatment); Prior stressors and no parasite exposure at 60 DPE
term estimate std.error statistic p.value p.adj.sig
Bray-Curtis_norm
(Intercept) −0.060 0.043 −1.398 0.165 ns
A+ T- P- 0.007 0.064 0.109 ≥0.25 ns
A- T+ P- 0.043 0.067 0.640 ≥0.25 ns
A+ T+ P- −0.119 0.079 −1.501 0.136 ns
Canberra_norm
(Intercept) −0.124 0.040 −3.122 0.002 **
A+ T- P- −0.089 0.059 −1.498 0.137 ns
A- T+ P- 0.099 0.062 1.602 0.112 ns
A+ T+ P- 0.123 0.073 1.680 0.096 ns
ANOVA of GLM - 60 DPE
ANOVA(GLM(Displacement ~ Treatment), type = 2); Prior stressors and no parasite exposure at 60 DPE
term statistic df p.value sig
Bray-Curtis_norm
Treatment 3.890 3.000 ≥0.25 ns
Canberra_norm
Treatment 12.030 3.000 0.007 **
Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); Prior stressors and no parasite exposure at 60 DPE
term .y. group1 group2 estimate std.error statistic adj.p.value Variable
Bray-Curtis_norm
Treatment Beta.Score A+ T- P- A- T- P- 0.007 0.064 0.109 ≥0.25 Treatment
Treatment Beta.Score A- T+ P- A- T- P- 0.043 0.067 0.640 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P- A- T- P- −0.119 0.079 −1.501 ≥0.25 Treatment
Treatment Beta.Score A- T+ P- A+ T- P- 0.036 0.069 0.514 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P- A+ T- P- −0.126 0.082 −1.545 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P- A- T+ P- −0.162 0.084 −1.931 0.213 Treatment
Canberra_norm
Treatment Beta.Score A+ T- P- A- T- P- −0.089 0.059 −1.498 ≥0.25 Treatment
Treatment Beta.Score A- T+ P- A- T- P- 0.099 0.062 1.602 ≥0.25 Treatment
Treatment Beta.Score A+ T+ P- A- T- P- 0.123 0.073 1.680 ≥0.25 Treatment
Treatment Beta.Score A- T+ P- A+ T- P- 0.187 0.064 2.921 0.018 Treatment
Treatment Beta.Score A+ T+ P- A+ T- P- 0.212 0.076 2.804 0.025 Treatment
Treatment Beta.Score A+ T+ P- A- T+ P- 0.025 0.078 0.321 ≥0.25 Treatment
Significant Results:Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); All treatments at 60 DPE
term .y. group1 group2 estimate std.error statistic adj.p.value Variable
Canberra_norm
Treatment Beta.Score A- T+ P- A+ T- P- 0.187 0.064 2.921 0.018 Treatment
Treatment Beta.Score A+ T+ P- A+ T- P- 0.212 0.076 2.804 0.025 Treatment
PCoA
Plot

Code
# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressNoParaExp_60DPE"



#### Dispersion --------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]] <-
  run_BetaDispersion(dist.matrix = beta.dist.mat[[tmp.resSubSection]], 
                     var = c("treatment_group"))

##### ANOVA ------------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA"]] <-
  get_HoD_anova(betaDisper = beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]],
                var = c("treatment_group"))

### Table

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA.Table"]] <-
  beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA"]] %>%
  
    # Create the table
    dplyr::group_by(Beta.Metric) %>%
    set_GT(var = "p.value", group.by = "Beta.Metric")  %>%
  
    # Title/caption
    gt::tab_header(
      title = "ANOVA: Homogeneity of Dispersion",
      subtitle = "ANOVA(Beta Disperson ~ Treatment); Prior stressors and no parasite exposure at 60 DPE"
    ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = c("p.value", "sig"),
      rows = p.value < 0.05
    )
  )

##### Tukey ------------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey"]] <-
  get_HoD_tukey(betaDisper = beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]],
                var = c("treatment_group")) 

### Table

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey.Table"]] <-
  beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey"]] %>%
    
    dplyr::mutate(
        group1 = case_when(
        group1 == "Control" ~ "A- T- P-",
        group1 == "Parasite" ~ "A- T- P+",
        group1 == "Antibiotics" ~ "A+ T- P-",
        group1 == "Antibiotics_Parasite" ~ "A+ T- P+",
        group1 == "Temperature" ~ "A- T+ P-",
        group1 == "Temperature_Parasite" ~ "A- T+ P+",
        group1 == "Antibiotics_Temperature" ~ "A+ T+ P-",
        group1 == "Antibiotics_Temperature_Parasite" ~ "A+ T+ P+",
    ),
        group2 = case_when(
        group2 == "Control" ~ "A- T- P-",
        group2 == "Parasite" ~ "A- T- P+",
        group2 == "Antibiotics" ~ "A+ T- P-",
        group2 == "Antibiotics_Parasite" ~ "A+ T- P+",
        group2 == "Temperature" ~ "A- T+ P-",
        group2 == "Temperature_Parasite" ~ "A- T+ P+",
        group2 == "Antibiotics_Temperature" ~ "A+ T+ P-",
        group2 == "Antibiotics_Temperature_Parasite" ~ "A+ T+ P+",
    )) %>%
  
  # Create the table
  dplyr::group_by(Beta.Metric) %>%
  set_GT(var = "adj.p.value", group.by = "Beta.Metric")  %>%
  
  # Title/caption
  gt::tab_header(
    title = "Tukey: Homogeneity of Dispersion",
    subtitle = "Tukey(Beta Disperson ~ Treatment); Prior stressors and no parasite exposure at 60 DPE"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = c("adj.p.value", "sig"),
      rows = adj.p.value < 0.05
    )
  )
Tables
ANOVA: Homogeneity of Dispersion
ANOVA(Beta Disperson ~ Treatment); Prior stressors and no parasite exposure at 60 DPE
term df sumsq meansq statistic p.value sig
bray
treatment_group 3.000 0.289 0.096 3.589 0.016 *
Residual 103.000 2.767 0.027 NA NA NA
canberra
treatment_group 3.000 0.076 0.025 6.573 <0.001 ***
Residual 103.000 0.395 0.004 NA NA NA
Tukey: Homogeneity of Dispersion
Tukey(Beta Disperson ~ Treatment); Prior stressors and no parasite exposure at 60 DPE
.y. term group1 group2 estimate conf.low conf.high adj.p.value sig
bray
Distance treatment_group A+ T- P- A- T- P- −0.036 −0.142 0.070 ≥0.25 ns
Distance treatment_group A- T+ P- A- T- P- −0.064 −0.174 0.046 ≥0.25 ns
Distance treatment_group A+ T+ P- A- T- P- −0.162 −0.294 −0.031 0.009 **
Distance treatment_group A- T+ P- A+ T- P- −0.028 −0.143 0.087 ≥0.25 ns
Distance treatment_group A+ T+ P- A+ T- P- −0.126 −0.262 0.009 0.077 ns
Distance treatment_group A+ T+ P- A- T+ P- −0.098 −0.237 0.040 ≥0.25 ns
canberra
Distance treatment_group A+ T- P- A- T- P- −0.056 −0.096 −0.016 0.002 **
Distance treatment_group A- T+ P- A- T- P- 0.010 −0.032 0.051 ≥0.25 ns
Distance treatment_group A+ T+ P- A- T- P- −0.023 −0.072 0.027 ≥0.25 ns
Distance treatment_group A- T+ P- A+ T- P- 0.066 0.023 0.109 <0.001 ***
Distance treatment_group A+ T+ P- A+ T- P- 0.034 −0.017 0.085 ≥0.25 ns
Distance treatment_group A+ T+ P- A- T+ P- −0.032 −0.085 0.020 ≥0.25 ns
CAP
Plot

Code
# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressNoParaExp_60DPE"


#### Capscale ----------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.mod"]] <-
  run_capscale(ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Treatment  %in% c(
          "A- T- P-",
          # "A- T- P+",
          "A+ T- P-",
          # "A+ T- P+",
          "A- T+ P-",
          # "A- T+ P+",
          "A+ T+ P-"#,
          # "A+ T+ P+"
          )) %>%
          microViz::ps_filter(Time == 60), 
               dist.matrix = beta.dist.mat[[tmp.resSubSection]], 
               formula_str = "dist ~ Treatment")

##### ADONIS ------------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS"]] <-
  run_cap_adonis(ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Treatment  %in% c(
          "A- T- P-",
          # "A- T- P+",
          "A+ T- P-",
          # "A+ T- P+",
          "A- T+ P-",
          # "A- T+ P+",
          "A+ T+ P-"#,
          # "A+ T+ P+"
          )) %>%
          microViz::ps_filter(Time == 60),
                 dist.matrix = beta.dist.mat[[tmp.resSubSection]], 
                 formula_str = "dist ~ Treatment",
                 by.method = "terms") 
## Warning in tidy.anova(.): The following column names in ANOVA output were not
## recognized or transformed: SumOfSqs, R2

## Warning in tidy.anova(.): The following column names in ANOVA output were not
## recognized or transformed: SumOfSqs, R2
### Table

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS.Table"]] <-
  beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS"]] %>%
    set_GT(var = "p.value", group.by = "Beta.Metric") %>%
  
  # Title/caption
  gt::tab_header(
    title = "ADONIS2",
    subtitle = "adonis2(Beta Distance ~ Treatment); Prior stressors and no parasite exposure at 60 DPE"
  )  %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = c("p.value", "sig"),
      rows = p.value < 0.05
    )
  )
Tables
ADONIS2
adonis2(Beta Distance ~ Treatment); Prior stressors and no parasite exposure at 60 DPE
term df SumOfSqs R2 statistic p.value sig
bray
Treatment 3.000 1.688 0.102 3.887 0.001 ***
Residual 103.000 14.910 0.898 NA NA NA
Total 106.000 16.598 1.000 NA NA NA
canberra
Treatment 3.000 3.157 0.086 3.234 0.001 ***
Residual 103.000 33.522 0.914 NA NA NA
Total 106.000 36.679 1.000 NA NA NA